Molecular hydrodynamics of the moving contact line in two-phase immiscible flows 
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The "no-slip" boundary condition, i.e., zero fluid velocity relative to the solid at the fluid-solid 
interface, has been very successful in describing many macroscopic flows, yet there is no compelling 
theoretical argument to justify this standard boundary condition of textbook continuum hydrody- 
namics. A problem of principle arises when the no-slip boundary condition is used to model the 
hydrodynamics of immiscible-fluid displacement in the vicinity of the moving contact line, where 
the interface separating two immiscible fluids intersects the solid wall. Decades ago it was already 
known that the moving contact line is incompatible with the no-slip boundary condition, since the 
latter would imply infinite dissipation due to a non-integrable singularity in the stress near the 
contact line. While subsequent molecular dynamics (MD) studies have clearly demonstrated fluid 
slipping relative to the wall at the contact line, the exact rule that governs this relative slip has 
eluded numerous prior attempts. In fact, over the years there have been numerous ad hoc mod- 
els and proposals aiming to resolve the incompatibility of the no-slip boundary condition with the 
moving contact line, but none was able to quantitatively account for the near-complete slip of the 
contact line observed in MD simulations. 

In this paper we first present an introductory review of the problem, including (1) the cause of the 
stress singularity, (2) the ad hoc introduction of the slip boundary condition, (3) the MD evidence of 
fiuid slipping, (4) the gap between the existing MD results and a continuum hydrodynamic descrip- 
tion, and (5) a preliminary account on how to bridge the MD results and a continuum description. 
We then present a detailed review of our recent results on the contact-line motion in immiscible 
two-phase fiow, from MD simulations to continuum hydrodynamics calculations. Through exten- 
sive MD studies and detailed analysis, we have uncovered the slip boundary condition governing 
the moving contact line, denoted the generalized Navier boundary condition. We have used this 
discovery to formulate a continuum hydrodynamic model whose predictions are in remarkable quan- 
titative agreement with the MD simulation results at the molecular level. These results serve to 
affirm the validity of the generalized Navier boundary condition, as well as to open up the possibility 
of continuum hydrodynamic calculations of immiscible fiows that are physically meaningful at the 
molecular level. 

PACS: 47.ll.-l-j, 68.08.-p, 83.10.Mj, 83.10.Ff, 83.50.Lh 



I. INTRODUCTION 

The no-slip boundary condition, i.e., zero relative tangential velocity between the fluid and solid at the interface, 
serves as a cornerstone in continuum hydrodynamics [1]. Although fluid slipping is generally detected in molecular 
dynamics (MD) simulations for microscopically small systems at high flow rate [2-5] , the no-slip boundary condition 
still works well for macroscopic flows at low flow rate. This is due to the Navier boundary condition which actually 
accounts for the slip observed in MD simulations [2-5]. This boundary condition, proposed by Navier in 1823 [6], 
introduces a slip length Ig and assumes that the amount of slip is proportional to the shear rate in the fluid at the 
solid surface: 

v^'^P = ~hn ■ [(Vv) + (Vv)^] , 
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where v''''^' is the slip velocity at the surface, measured relative to the (moving) wall, [(Vv) + (Vv)-^] is the tensor 
of the rate of strain, and n denotes the outward surface normal (directed out of the fluid). According to the Navier 
boundary condition, the slip length is defined as the distance from the fluid-solid interface to where the linearly 
extrapolated tangential velocity vanishes (see Fig. 1). Typically, Ig ranges from a few angstroms to a few nanometers 
[2-5]. For a flow of characteristic length R and velocity U, the slip velocity is on the order of Ulg/R. This explains why 
the Navier boundary condition is practically indistinguishable from the no-slip boundary condition in macroscopic 
flows where v^^'p /U ~ Is/R 0. 




solid 

FIG. 1. Slip length introduced in the Navier boundary condition, defined as the distance from the fluid-solid interface to 
where the linearly extrapolated tangential velocity vanishes. 

It has been well known that the no-slip boundary condition is not applicable to the moving contact line (MCL) 
where the fluid-fluid interface intersects the solid wall [7-9] (see Fig. 2 for both the static and moving contact lines). 
The problem may be simply stated as follows. In the two-phase immiscible flow where one fluid displaces another 
fluid, the contact line appears to "slip" at the solid surface, in direct violation of the no-slip boundary condition. 
Furthermore, the viscous stress diverges at the contact line if the no-slip boundary condition is applied everywhere 
along the solid wall. This stress divergence is best illustrated in the reference frame where the fluid-fluid interface 
is time-independent while the wall moves with velocity U (see Fig. 2b). As the fluid velocity has to change from U 
at the wall (as required by the no-slip boundary condition) to zero at the fluid-fluid interface (which is static), the 
viscous stress varies as rjU /x^ where rj is the viscosity and x is the distance along the wall away from the contact line. 
Obviously, this stress diverges as a; ^ because the distance over which the fluid velocity changes from U to zero 
tends to vanish as the contact line is approached. In particular, this stress divergence is non-integrable (the integral 
of 1/a; yields In a;), thus implying infinite viscous dissipation. 
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FIG. 2. (a) Static contact line. A fluid-fluid interface is formed between two immiscible fluids and intersects the solid wall 
at a contact line. The static contact angle Os is determined by the Young's equation: 7 cos 0s -I- 72 — 71 = 0, where 71, 72, and 
7 are the three interfacial tensions at the three interfaces (two fluid-solid and one fluid-fluid), (b) Moving contact line. When 
one fluid displaces another immiscible fluid, the contact line is moving relative to the solid wall. (Here fluid 1 displaces fluid 
2.) Due to the contact-line motion, the dynamic contact angle dd deviates from the static contact angle ds- We will show that 
this deviation is primarily responsible for the near-complete slip at the contact line. 
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Over the years there have been numerous models and proposals aiming to resolve the incompatibility of the no-slip 
boundary condition with the MCL. For example, there have been the kinetic adsorption/desorption model by Blake 
[10], the fluid slip models by Hocking [11], by Huh and Mason [12], and by Zhou and Sheng [13], and the Cahn-Hilliard- 
van der Waals models by Seppecher [14], by Jacqmin [15], and by Chen et al. [16]. In the kinetic adsorption/desorption 
model by Blake [10], the role of molecular processes was investigated. A deviation of the dynamic contact angle from 
the static contact angle was shown to be responsible for the fluid/fluid displacement at the MCL. In the slip model 
by Hocking [11], the effect of a slip coefficient (slip length) on the flow in the neighborhood of the contact line was 
examined. Two slip models were used by Huh and Mason [12]: Model I for classical slippage assumes the slip velocity 
is proportional to the shear stress exerted on the solid; Model II for local slippage assumes that over a (small) distance 
the liquid slips freely where fluid stress vanishes, but thereafter the liquid/solid bonding has been completed and the 
no-slip boundary condition is applied. In the slip model by Zhou and Sheng [13], the incompressible Navier-Stokes 
equation was solved using a prescribed tangential velocity profile as the boundary condition, which exponentially 
interpolates between the complete-slip at the contact line and the no-slip far from the contact line. The Cahn- 
Hilliard-van der Waals models by Seppecher [14], by Jacqmin [15], and by Chen et al. [16] suggested a resolution 
when perfect no-slip is retained. With the fluid-fluid interface modeled to be diffuse, the contact line can thus move 
relative to the solid wall through diffusion rather than convection. All the above models are at least mathematically 
valid because the divergence of stress has been avoided, citlica' by introducing some molecular process to drive the 
contact line [10], or by allowing slip to occur [11-13], or by modeling a diffuse fluid-fluid interface [14-16]. Pismen 
and Pomeau have presented a rational analysis of the hydrodynamic phase field (diffuse interface) model based on 
the lubrication approximation [17]. 

The most usual (and natural) approach to resolve the stress divergence has been to allow slip at the solid wall 
close to the contact line. In fact, molecular dynamics (MD) studies have clearly demonstrated the existence of fluid 
slipping in the molecular-scale vicinity of the MCL [18,19]. However, the exact rule that governs this relative slip 
has eluded numerous prior attempts. As a matter of fact, none of the existing models has proved successful by 
quantitatively accounting for the contact-line slip velocity profile observed in MD simulations. In a hybrid approach 
by Hadjiconstantinou [20], the MD slip profile was used as the boundary condition for finite-element continuum 
calculation. The continuum results so obtained match the corresponding MD results, therefore demonstrating the 
feasibility of hybrid algorithm [21,22]. But the problem concerning the boundary condition governing the contact-line 
motion was still left unsolved. 

In this paper we first present an introductory review of the problem, including (1) the origin of the stress singularity, 
(2) the ad hoc introduction of the slip boundary condition, (3) the MD evidence of fluid slipping, (4) the gap between 
the existing MD results and a continuum hydrodynamic description, and (5) a preliminary account on how to bridge 
the MD results and a continuum description. We then present a detailed review of our recent results on the contact- 
line motion in immiscible two-phase flow, from MD simulations to continuum hydrodynamics calculations [23]. In 
our MD simulations, we consider two immiscible dense fluids of identical density and viscosity, with the temperature 
controlled above the gas-liquid critical point. (Similar results would be obtained if the temperature is below the 
critical point.) For fluid-solid interactions, we choose the wall density and interaction parameters to make sure (1) 
there is no epitaxial locking of fluid layer(s) to the solid wall, (2) a finite amount of slip is allowed in the single- 
phase flow for each of the two inuriiscible fluids, (3) the fluid-fluid interface makes a finite microscopic contact angle 
with the solid surface. Through extensive MD simulations and detailed analysis, we have uncovered the boundary 
condition governing the fluid slipping in the presence of a MCL. With the help of this discovery, we have formulated 
a continuum hydrodynamic model of two-phase immiscible flow. Numerical solutions have been obtained through an 
explicit flnite-difference scheme. A comparison of the MD and continuum results shows that velocity and fluid-fluid 
interface profiles from the MD simulations can be quantitatively reproduced by the continuum model. 

The paper is organized as follows. In Sec. II we review a few known facts in mathematics and physics concerning 
the contact line motion. Together, they point out the right direction to approach and elucidate the problem. In 
fact, they almost tell us what is expected for a boundary condition governing the slip at the MCL and in its vicinity. 
In Sec. HI we outline the main results from MD simulations to continuum hydrodynamic modeling. In Sec. IV 
we present the first part of the MD results. From the slip velocity and the tangential wall force measured at the 
fluid-solid interface, a slip boundary condition is deduced. In Sec. V we formiilate a continuum hydrodynamic model 
of two-phase immiscible flow. The continuum differential expression for the tangential stress at the solid surface is 
derived, from which the generalized Navier boundary condition (CNBC) is obtained from the slip boundary condition 
deduced in Sec. IV. In Sec. VI wo show a systematic comparison of the MD and continuum results. The validity 
of the GNBC and the continuum model is demonstrated. In Sec. VII we present the second part of the MD results, 
concerning the tangential force balance in a boundary layer at the fluid-solid interface and the decomposition of the 
tangential stress in the fluid-fluid interfacial region. In Sec. VIII we establish the correspondence between the stress 
components measured in MD and those defined in the continuum hydrodynamics. Based on this correspondence, the 
continuum GNBC is obtained in an integrated form from the MD results in Sees. IV and VII. This may be regarded 
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as a direct MD verification of the continuum GNBC. It also justifies tire use of the Cahn-Hilliard hydrodynamic 
formulation of two-phase flow, from which the continuum form of GNBC, as verified by the MD results, is derived. 
The paper is concluded in Sec. IX. 

II. STRESS AND SLIP: A BRIEF REVIEW 

A. Non-integrable stress singularity: the Huh-Scriven model 

Hua and Scriven [7] considered a simple model for the immiscible two-phase flow near a MCL. A flat solid surface is 
moving with steady velocity U in its own plane and a flat fluid-fluid interface, formed between two immiscible phases 
A and B, intersects the solid surface at a contact line (see Fig. 3). The contact line is taken as the origin of a polar 
coordinate system {r,9), in which the contact angle is <j). 




U 

FIG. 3. The Hua-Scriven model. 



The two-dimensional flow of Newtonian and incompressible fluids is dominated by viscous stress for r ^ v/pU, 
where r] is the viscosity and p the mass density. In the viscous flow approximation, the Navicr-Stokes equation is 
linearized and the steady flow is solved from the biharmonic equation for the stream function tp{7\ 9): 

v2(vV) = 0. 

The boundary conditions to be imposed at = 0, tt (solid surface) and 9 = (p (fluid-fluid interface) are: (i) vanishing 
normal component of velocity at the solid surface and fluid-fluid interface, (ii) continuity of tangential velocity at 
the fluid-fluid interface, (iii) continuity of tangential viscous stress at the fluid-fluid interface, (iv) no slip at the 
solid surface. Here conditions (i)-(iii) are well justifled while (iv) is usually taken as a postulate of continuum 
hydrodynamics. The no-slip boundary condition can be justified a posteriori in macroscopic fiow by checking the 
correctness of its consequences. In the present model, however, it leads to physically incorrect results for stress, in 
the microscopic vicinity of the MCL. 

The similarity solution of the biharmonic equation is in the form of 

ip{r, 6) = r{asm6 + bcosO + c^sin^ -|- d6 cos 6), 

in which the eight constants (4 for phase A and 4 for phase B) are determined by the eight boundary conditions in 

(i)-(iv). What Hua and Scriven found is that the shear stress and pressure fields vary as l/r and hence increase in 
magnitude without limit as the contact line r = is approached. As a consequence, the total tangential force exerted 
on the solid surface, which is an integral of the tangential stress along the surface, is logarithmically infinite. That 
indicates a singularity in viscous dissipation, which is physically unacceptable. 

Obviously, the Hua-Scrivcn model is defective in the immediate vicinity of the MCL. This is also seen from the 
normal stress difference across the fluid-fluid interface, which varies as l/r as well. According to the Laplace's 
equation, the interface curvature should increase rapidly as the contact line (r = 0) is approached, in order to balance 
the normal stress difference by curvature force. This is clearly inconsistent with the assumption of a flat fluid-fluid 
interface. Nevertheless, the flow field solved from the Hua-Scriven model may approximately describe the asymptotic 
region at a large distance from the contact line (where the viscous flow approximation is still valid). In that region, 
the no-slip boundary condition is considered valid and the fluid-fluid interface is almost flat, due to the reduced stress. 

B. Introducing the slip boundEiry condition 

The deflciency of the Hua-Scriven model results from the incompatibility of the no-slip boundary condition with a 
MCL: no slip means Vr = ±U a.t the solid surface {6 = 0, n) where r > while at r = the MCL requires a perfect 
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slip. That is, the no-sHp boundary condition leads to a velocity discontinuity at the MCL. In order to remove the 
stress singularity at the MCL, while retaining the Newtonian behavior of stress, the continuity of velocity field must 
be restored. For this purpose, a slip profile can be introduced to continuously interpolate between the complete slip 
at the MCL and the no-slip boundary condition that must hold at regions far away. 

Dussan V. [24] considered a plate of infinite extent either inserted into or withdrawn from a semi-infinite domain 
of fluid at a constant velocity Uq (see Fig. 4). The contact line is taken as the origin of a polar coordinate system 
(r, (j)), in which the apparent contact angle is a at r — > oo. The equation of motion is the Navier-Stokes equation with 
the incompressibility condition V • u = 0. The boundary conditions at the solid surface (p = and the free surface 
{(r, (/!)(r)}|0 < r < oo} are: (i) the kinematic boundary conditions u-^ = Oat(/) = and u • n = at = (t>(r), where 
n is an outward unit vector normal to the free surface, (ii) the dynamic boundary condition t-T-n = Oat^ = ^(r), 
where T is the stress tensor and t is a unit vector tangent to the free surface, i.e., t • n = 0, (iii) the Laplace condition 
n • T • n = an at — (t){r), where a is the surface tension and n the interface curvature, (iv) the slip boundary 
condition u = U{r)r at = 0, where U{r) is a prescribed function, and (v) the critical static contact angle (/>(0) = (j)s- 



Free surface 




6=0 



FIG. 4. A plate is inserted into a semi-infinite domain of fluid at a constant velocity Uo- The angle between the plate and 
the free surface at r — > oo is a. 



The slip boundary condition must continuously interpolate between the complete slip at the MCL (r 
no-slip boundary condition at r ^ oo: 



0) and the 



lim U{r) = 0, 

r— »0 



lim U{r) = Uq. 



However, the form of U{r) in the intermediate region, where U varies from to Uq, is unknown. Three different slip 
boundary conditions were used for U{r), in order to assess the sensitivity of the overall flow field to the form of the 
slip boundary condition. They are 



Ui 



r/L, 



1 + r/R 



■Uo, U2 



1 + {r/L,y 



:Uo, Ui 



l + (r/L.)V2^°- 



where is a length scale called the slip length (not the one defined in the Navier boundary condition) . It was found 
that on the slip length scale the fiow fields are quite different whereas on the meniscus length scale, i.e., the length 
scale on which almost all fluid-mechanical measurements are performed, all the flow fields are virtually the same. 
That is, identical macroscopic flow behaviors are expected from different slipping models. 



C. Slip observed in molecular dynamics simulations 

According to the conclusion in Ref. [24], only in the microscopic slip region (of length scale ^ Ls) can different slip 
models be distinguished. This makes the experimental verification of a particular slip model very difficult because 
experiments usually probe distances (~ 1/xm) much larger than Lg- Naturally, computer experiments become very 
useful in elucidating the problem [25]. 

Non-equilibrium MD simulations were carried out to investigate the fluid motion in the vicinity of the MCL, in both 
the Poiseuille-flow [18] and Couette-flow [19] geometries. In these MD simulations, interactions between fluid molecules 
were modeled with the Lcnnard- Jones potential, modified to segregate immiscible fluids. The confining walls were 
constructed with a molecular structure. Wall-fluid interactions were also modeled with the Lennard-Jones potential, 
with energy and length scales different from those of the fluid-fluid interactions. In the simulations performed in the 
Couettc geometry [19], two immiscible fluids were confined between two planar walls parallel to the xy plane and a 
shear fiow was induced by moving the two solid walls in ±x directions at constant speed U (see Fig. 5). Steady-state 
velocity fields were obtained from the time average of fluid molecular velocities in small bins. 
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FIG. 5. The xz projection of the simulated system. Thick soUd Unes represent soUd walls and dashed lines represent fluid-fluid 
interfaces. The arrows indicate fluid velocity close to the solid walls, from which the variation of the amount of slip is clearly 
seen. 



There was clean evidence for a slip region in the vicinity of the MCL: appreciable slip occurs within a length scale 
~ lOcr, where a is the length scale in the fluid-fluid Lennard-Jones potential, and at the MCL the slip is near-complete. 

Far from the interfaces, viscous damping makes the flow insensitive to the fast variations near the interfaces, and 
hence uniform shear flow was observed, with a negligible amount of slip. Therefore, MD simulations provide evidence 
for a cutoff, below which the no-slip boundary condition breaks down, introduced phenomenologically in various slip 
models to remove the stress singularity. However, the exact boundary condition that governs the observed slip was 
left unresolved. In particular, a breakdown of local hydrodynamics in the molecular scale slip region was suggested 
[19], considering the extreme velocity variations there. 

The Navier slip model assumes that the amount of sUp is proportional to the tangential component of the stress 
tensor, P^zj in the fluid at the solid surface. In Ref. [19], the microscopic value of Pxz was directly measured. A 
comparison to the slip profile is roughly consistent with the Navier slip model. However, a large discrepancy was 
found between the microscopic value of P^z and the shear rate dv^/dz. According to the authors, this discrepancy 
suggests a breakdown of local hydrodynamics. On the contrary, we will show in this paper: 

(i) The Navier slip model is the correct model describing the fluid slipping near the MCL. 

(ii) The tangential stress in the Navier model is not merely viscous. 

(iii) The interfacial tension plays an important role in governing the contact-line slip in immiscible two-phase flows. 

(iv) There is no breakdown of local hydrodynamics. 



D. Fluid/fluid displacement driven by unbalanced Young force 

1. Unbalanced Young force 

For a cylindrical capillary of radius r, if the steady displacement is sufficiently slow, then the pressure drop across the 
moving fluid-fluid interface is given by Ap = 2712 cos 0/r, where 712 is the interfacial tension, and 9 is an appropriate 
dynamic contact angle. At equilibrium, the pressure drop is given by Ap° = 2712 cosO^/r, where 9° is the equilibrium 
contact angle. In general, 9 and 9o differ. Physically, Ap'^ measures the change of the interfacial free energy at the 
fluid-solid interface when the fluid-fluid interface moves relative to the wall, and Ap measures the external work done 
to the system. Therefore, the difference Ap — Ap^ is a measure of the dissipation due to the presence of the MCL. 

Let's consider a displacement of fluid 2 by fluid 1 over distance L (see Fig. 6). According to the Young's equation 
for the equilibrium contact angle, the force Trr^Ap*^ equals 27rr7i2 cosfl*^ ~ 2Trr{'jsi — 7S2), where 751 and 752 are the 
interfacial tensions for the 5/1 and S/2 interfaces, respectively. Thus the change of the interfacial free energy at the 
fluid-solid interface is given by wr'^Ap^L = 27rri(75i — 752)- The external work done to the system is simply nr^ApL. 
It follows that the dissipation due to the presence of the MCL is given by 7rr^(Ap — Ap^)L = 27rrL7i2 (cos ^ — cos 6°), 
where 712 (cos 9 — cos 9'^) — 712 cos 9 + 752 — 7si is the unbalanced Young force [9]. 
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FIG. 6. A displacement of fluid 2 by fluid 1 over distance L in a cylindrical capillary of radius r. 
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In order to find a relation between the displacement velocity v and the unbalanced Young force Fy, two classes 
of models have been proposed to describe the contact line motion: a) An Eyring approach, based on the molecular 
adsorption/desorption processes at the contact line [10]. b) A hydrodynamic approach, assuming that dissipation is 
dominated by viscous shear flow inside the wedge [9] . Viscous flows in wedges were calculated by Hua and Scriven [7] . 
For wedges of small (apparent) contact angle, a lubrication approximation can be used to simplify the calculations 
[9]. As discussed in Sees. II A and II B, a (molecular scale) cutoff has to be introduced to remove the logarithmic 
singularity in viscous dissipation. On the contrary, the Eyring approach assumes that the molecular dissipation at 
the tip is dominant. A quantitative theory is briefly reviewed below. 



2. Blake's kinetic model 



The role of interfacial tension was investigated in a kinetic model by Blake and Haynes [10] . Consider a fluid- fluid 
interface in contact with a flat solid surface at a line of three-phase contact (see Fig. 7a). Viewed on a molecular 
scale, the three-phase line is actually a fluctuating three-phase zone, where adsorbed molecules of one fluid (at the 
solid surface) interchange with those of the other, either by migration at the solid surface or through the contiguous 
bulk phases. In equilibrium the net rate of exchange will be zero. 

For a three-phase zone moving relative to the solid wall (Fig. 7a), the net displacement, driven by the unbalanced 
Young force Fy = 712 cosd + 752 — 751, is due to a nonzero net rate of exchange. Let ^ be the interfacial thickness, 
a be the area of an adsorption site, and A be the hopping distance of molecules. The force per adsorption site is 
approximately crFy/$, and the energy shift over distance A is approximately XaFy/S, ~ FyX"^ (see Fig. 7b). This 
energy shift leads to two different rates and K- : 



K± = k exp 



for forward and backward hopping events, respectively. Here W is an activation energy for hopping. It follows that 
the velocity of the MCL is v = X{K+ - K_) oc sinh [FyX^ /2kBT) . For very small Fy,v(xFy. 




FIG. 7. (a) A molecular picture of the three-phase zone, (b) Shifted potential profile. 

Blake's kinetic model shows that fluid slipping can be induced by the unbalanced Young force at the contact line. 
Therefore, it emphasizes the role of interfacial tension, though not in a hydrodynamic formulation. On the contrary, 
the authors of Ref. [19] considered the viscous shear stress as the only driving force. In fact a large discrepancy was 
found between the shear rate and the microscopic value of the tangential stress (the driving force according to the 
Navier slip model). In this paper we will show that in the two-phase interfacial region, such a discrepancy is exactly 
caused by the neglect of the non- viscous contribution from interfacial tension. 

E. From the Navier boundciry condition to the generalized Navier boundary condition: a preliminary 

discussion 

Here we give a preliminary account on the main finding reported in Ref. [23], the GNBC. Based on the results 
reviewed in Sees. II A, II B, II C, and II D, we try to show: a) The Navier boundary condition may actually work for 
the single-phase slip region near the MCL. b) In the two-phase contact-line region, the GNBC is a natural extension 
of the Navier boundary condition, with the fluid-fluid interfacial tension taken into account. 
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1. Navier boundary condition and slip length 

The validity of the Navier boundary condition v*'*^ = —Ign ■ [(Vv) + (Vv)^] has been well established by many 
MD studies for single-phase fluids [2 5]. This boundary condition is a constitutive equation that locally relates the 
amount of slip to the shear rate at the solid surface, though in most of the reported simulations, the hydrodynamics 
involves no spatial variation along the solid surface. Physically, the Navier boundary condition is local in nature 
simply because intermolecular wall-fluid interactions are short-ranged. 

The slip length Is is a phenomenological parameter that measures the local viscous coupling between fluid and solid. 
Figure 8 illustrates the viscous momentum transport between fluid and solid through a boundary layer of fluid. The 
thickness of this boundary layer, zq, must be of molecular scale, within the range of wall-fluid interactions. Now we 
show that the slip length Is is defined based on a linear law for tangential wall force and the Newton's law for shear 
stress. 



FIG. 8. A boundary layer of fluid at the fluid-solid interface, responsible for viscous momentum transport between fluid and 
solid. Circles represent fluid molecules and solid squares represent wall molecules. 

Hydrodynamic motion of fluid at the solid surface is measured by the slip velocity u*'*'' in the x direction, defined 
relative to the (moving) wall. When is present, a tangential wall force is exerted on the boundary- layer 
fluid, defined as the rate of tangential {x) momentum transport per unit wall area, from the wall to the (boundary- 
layer) fluid. Physically, this force represents a time-averaged eflfect of wall-fluid interactions, to be incorporated into a 
hydrodynamic slip boundary condition. The linear law for is expressed by G!^ = ~(3v^''^^, where f3 is called the slip 
coefficient and the minus sign means the fluid-solid c;oupling is viscous (frictional). For the boundary layer of molecular 
thickness, the tangential wall force must be balanced by the tangential fluid force Gl = dz {dx(Jxx + dzCTzx), 
where the integral is over the z-span of the boundary layer (see Fig. 8), a^xizx) are the xx{zx) components of the 
fluid stress tensor, and dx<Jxx + dz<^zx is the fluid force density in the x direction. Prom the equation of tangential 
force balance -|- G^ = 0, we obtain 



I dz {dxCTxx + dz(Tzx) = dx dza^xix, z) + a^xix, zq) = /3v*''P(a;). 

Jo 

This equation should be regarded as a microscopic expression for the Navier slip model: the amount of slip is 
proportional to the tangential fluid force at the fluid-solid interface. When the normal stress axx varies slowly in the 
X direction and the tangential stress a^x is caused by shear viscosity only, the above equation becomes rjdzVxix, zq) = 
Pv^'-^P{x), where r] is the viscosity and dzVx{x,zo) is the shear rate "at the solid surface". For flow fields whose 
characteristic length scale is much larger than zo, the boundary layer thickness, we replace zq hy z = and recover 
the Navier boundary condition for contimmm hydrodynamics, in which the slip length is given by Ig = rj/p. To 
summarize, the hydrodynamic viscous coupling between fluid and solid is actually measured by the slip coefficient /3. 
The Navier boundary condition, in which the slip length is introduced, is due to the Newton's law for viscous shear 
stress. For very weak viscous coupling between fluid and solid, /3 ^ 0, and thus Ig — > oo, making dzVx 0: the 
fluid-solid interface becomes a free surface. 



2. Single-phase region 



A number of MD studies have shown that for single-phase flows, the Navier boundary condition is valid in describing 
the fluid slipping at solid surface [2-5]. Therefore, we expect that it can also describe the partial slip observed in the 
single-phase region near the MCL [19]. However, according to the authors of Ref. [19], the Navier boimdary condition 
failed even in the single-phase region: the velocity gradient dvx/dz was not proportional to the amount of slip. This 
discrepancy was attributed to a breakdown of local hydrodynamics, considering the very large velocity variations 
observed near the MCL. Here we present a heuristic discussion, to explain why such a discrepancy is expected even 
if the Navier boundary condition actually works for the single-phase region. 
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First, the success of the hybrid approaches outUncd below strongly indicates that local hydrodynamics doesn't break 
down in the slip region. In Ref. [19], an apparent contact angle Oapp was defined at half the distance between two 
solid surfaces (see Fig. 5). For Oapp < 135°, fluid-fluid interfaces could be approximated by planes. The Navier-Stokes 
equation was solved for the simplified geometry. For this purpose, the tangential velocity along fluid-fluid interface 
was set to be zero (according to the simulation results) and the slip velocity Ay(r) (measured relative to the moving 
wall) was specified as a function of distance r from the MCL: Ay(r) = ±[/exp[— rln2/S'], with -I- and — for the 
lower and upper walls, respectively. This Ay(r), proposed according to the MD slip profile, continuously interpolates 
between the complete slip at the MCL (AF(r) = ±XJ at r = 0) and the no-slip boundary condition {/S.V{r) ^ at 
r oo). With a proper value chosen for S (~ l.Sfx), the MD flow fields were reproduced by solving the Navier-Stokes 
equation with the above boundary conditions. Recently, an improved hybrid approach has been used to reproduce 
the MD simulation results for contact-line motion in a Poiseuille geometry [20]. 

To further test the validity of continuum modeling, we have solved the Navier-Stokes equation for a corner flow 
[26]. Consider one rigid plane sliding steadily over another, with constant inclination Oapp (see Fig. 9). The fluid 
is Newtonian and incompressible. The no-slip boundary condition is applied at the vertical plane and the Navier 
boundary condition applied at the horizontal plane. The kinematic condition of vanishing normal component of 
velocity at the solid surface requires = at the vertical plane and t^z = at the horizontal plane, and hence v = 
at the intersection O, which is taken as the origin of the coordinate system {x, z). This corner-flow model resembles 
the continuum model used in the hybrid approach above, and produces similar flow fields for quantitative comparison 
with MD results. In particular, the corner flow exhibits a slip profile very close to AF(r): the slip velocity + U 
shows a linear decrease over a length scale ^ Ig, the slip length in the Navier condition, followed by a more gradual 
decrease. (Note that + U = U at O implies complete slip.) 



kZ 



e,pp=jt/2 



no-sIip condition 

Navier-Stokes equation 

Navier condition 



FIG. 9. Two-dimensional corner flow caused by one rigid plane sliding over another, with constant inclination Oapp = 90°. 
In the reference frame of the vertical plane, the horizontal plane is moving to the left with constant speed U. 



From the hybrid approach with the prescribed slip velocity AV{r) to the corner-flow model with the Navier boundary 
condition, they indicate that the single-phase region near the MCL can be modeled by the Navier-Stokes equation for 
an incompressible Newtonian fluid, supplemented by appropriate boundary conditions. Then a simple question arises: 
Given a continuum hydrodynamic model that uses the Navier boundary condition and approximately reproduces the 
MD flow fields in the single-phase region near the MCL, why do the MD simulation results show a large discrepancy 
between the velocity gradient dv^/dz and the amount of slip? 

The answer lies in the fast velocity variation in the vicinity of the MCL, where the flow is dominated by viscous 
stress. With the characteristic velocity scale set by U and the characteristic length scale set by the normal stress (Jxx 
must show as well a fast variation along the solid surface: dxi^xx ~ vU/l'j.. According to Sec. HE 1, the microscopic 
value of the tangential fluid force is given by dx J^" dzaxxix, z) + Ozxix, zo). This expression is necessary because 
G| is distributed in a boundary layer of molecular thickness zq. Obviously, to represent by (j^x ~ rjdvx/dz at 
z = zq alone, the normal stress cr^x near the solid surface has to vary slowly in the x direction. This is not the case 
in the vicinity of the MCL if the Navier slip model works there, as implied by the continuum corner-flow calculation 
which yields dxCFxx ^ vU/lg near the intersection. (It is reasonable to expect that if the Navier slip model is valid, then 
the normal stress measured in MD simulations should in general agree with that from continuum model calculation 
with the Navier condition.) Given dx Jq" dzaxx ~ vUzq/II according to the corner-flow model, and that zq and I a 
are both ^ la, it is obvious that considering only rjdvx/dz at z = zq would lead to an appreciable underestimate of 
Gl ~ rjU/ls in the slip region. In short, the large discrepancy between the velocity gradient dvx/dz and the amount 
of slip, observed in the single- phase slip region near the MCL, cannot be used to exclude the microscopic Navier slip 
model and the hydrodynamic Navier slip boundary condition, for if the Navier model is valid, then the tangential 
viscous stress rjdvx/dz, measured at some level away from the solid surface, is not enough for a complete evaluation 
of the tangential fluid force G^. 
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3. Two-phase region 



The linear law for tangential wall force, = — describes the hydro dynamic viscous coupling between 
fluid and solid. Assume that in the two-phase region, the two fluids interact with the wall independently. Then 
the tangential wall force becomes G'^ = —j^v^^^ at the contact line, with the slip coefficient given by /3 = {(3iPi + 
l3iP2)/{pi + P2), where (5\ and /32 are the slip coefBcients for the two single-phase regions separated by the fluid- fluid 
interface, and pi and p2 are the local densities of the two fluids in the contact-line region. Obviously, /3 varies between 
(3i and across the fluid-fluid interface. 

The equation of tangential force balance -\- G^. = Q must hold as well in the two-phase region. Therefore the 
Navier slip model is still of the form 

dx / dzaxx{x, z) -\- azx{x, zq) = /3v*''^(a;). 

Nevertheless, this does not lead to the Navier boundary condition rjdzVx = fiv^^^P anymore because in the two-phase 
region, the tangential stress is not contributed by shear viscosity only: there is a non- viscous component in Ozx, 
caused by the fluid-fluid interfacial tension. Put in an ideal form of decomposition, the tangential stress azx{x-,z) at 
level z can be expressed as 

(^zx{x, z) = ctI^{x, z) + a^^{x, z), 

where a^x is the viscous component due to shear viscosity and cr^ (the tangential Young stress) is the non-viscous 
component, which is narrowly distributed in the two-phase interfacial region and related to the interfacial tension 
through J._^^dxaJx{x,z) — j cos 9{z). Here J^^fdx denotes the integration across the fluid-fluid interface along the 
X direction, 7 is the interfacial tension, and 9{z) is the angle between the interface and the x direction at level z. 
Physically, the existence of a fluid-fluid interface causes an anisotropy in the stress tensor in the two-phase interfacial 
region. The interfacial tension is an integrated measure of that stress anisotropy. When expressed in a coordinate 
system different from the principal system, the stress tensor is not diagonalizcd and a nonzero cr^ appears. In the 
presence of shear flow, while the shear viscosity leads to the viscous component in a^x ) the non- viscous component 
is still present. A detailed discussion on the stress decomposition will be given in Sec. VII C. 
Consider an equilibrium state in which the tangential stress a^x is balanced by the gradient of the normal stress 
a° : 

^ XX 

dx / dza^xix, z) + a-1x{x, zq) = 0, 
Jo 

where the superscript denotes equilibrium quantities. Here the equilibrium tangential stress cr^^, is narrowly dis- 
tributed in the two-phase interfacial region and related to the interfacial tension through /.^^ dxcj^^ {x-i -z) = 7 cos 6^{z). 

(There is no viscous stress in equilibrium, and hence (t°^ is caused by the interfacial tension only.) Subtracting the 
equation of equilibrium force balance from the expression of Navier slip model, we obtain 

dxio^ dz[axx{x,z) - cr°^(a;,z)] + [a^^(a;,zo) - (J°x(.x,Zq)] 
= dx dz[axx{x, z) - alx{x, z)] + [al^ix, zq) + aj^ix, Zq) - a°x{x, Zq)] 
= l3v^^'P{x). 

This equation will be the focus of our continuum deduction from molecular hydrodynamics. In fact it leads to the 
GNBC which governs the fluid slipping everywhere, from the two-phase contact-line region to the single-phase regions 
away from the MCL. This will be elaborated in Sees. IV, V, VII, and VIII. 

To summarize, our preliminary analysis shows that compared to the single-phase region where the tangential stress 

is due to shear viscosity only, the two- phase region has the tangential Young stress as the additional component. 
Naturally, the Navier boundary condition, which considers the tangential viscous stress only, needs to be generalized 
to include the tangential Young stress. 

III. STATEMENT OF RESULTS 

We have carried out MD simulations for immiscible two-phase flows in a Couette geometry (see Fig. 10) [23]. The 
two immiscible fluids were modeled by using the Lennard-Jones (LJ) potentials for the interactions between fluid 
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molecules. The solid walls were modeled by crystalline plates. More technical details will be presented in Sees. IV 
and VII. The purpose of carrying out MD simulations is threefold: (1) To uncover the boundary condition governing 
the MCL, denoted the GNBC; (2) To fix the material parameters (e.g. viscosity, interfacial tension, etc) in our 
hydrodynamic model; (3) To produce velocity and interface profiles for comparison with the continuum hydrodynamic 
solutions. 



V- 



61 ■ ■ ■ ■■ ■ ■■ ■■■■■■ m,m 
o^o Oq qqOqO oq!; 

Bo8o°§o°ooo. 
o Qpo _o oO_o^^ -• •'^ •S \ 

Pb bY b b f b b III I 



O qO 0_0r~ 



I I I • I I U I J 



III ri* 



V 



FIG. 10. Scliematic of tlie immiscible Couette flow. Two immiscible fluids (empty and solid circles) are confined between 
two solid walls (solid squares) that are parallel to the xy plane and separated by a distance H. The Couette flow is generated 
by moving the top and bottom walls at a speed V along ±x, respectively. 



Our main finding is the GNBC: 



_ ~ _ „v , -Y 



Here (3 is the slip coefficient, v^'^^ is the tangential slip velocity measured relative to the (moving) wall, and azx 
is the hydrodynamic tangential stress, given by the sum of the viscous stress cr"^ and the uncompensated Young 
stress a^. The validity of the GNBC has been verified by a detailed analysis of MD data (Sees. IV, VII, and VIII) 
plus a comparison of the MD and continuum results (Sec. VI). Compared to the conventional Navier boundary 
condition which includes the tangential viscous stress only, the GNBC captures the uncompensated Young stress 
as the additional component. Together, the tangential viscous stress and the uncompensated Young stress give rise 
to the near-complete slip at the MCL. The uncompensated Young stress arises from the deviation of the fiuid-fluid 
interface from its static configuration and is narrowly distributed in the fluid-fluid interfacial region. Obviously, far 
away from the MCL, the uncompensated Young stress vanishes and the GNBC becomes the usual Navier Boundary 
condition (3vf'^^ = a^^ . 

We have incorporated the GNBC into the Cahn-Hilliard (CH) hydrodynamic formulation of two-phase flow [15,16] 
to obtain a continuum hydrodynamic model [23]. The continuum model may be briefly described as follows. The 

CH free energy functional [27] is of the form i^[(/)] = J dr K {'V(f>)'^ /2 + /(^) , where (j>{r) is the composition field 

defined by ^(r) = {p2 — pi)/{p2 + Pi)j with pi and p2 being the local number densities of the two fluid species, 
/(0) = —r(jP'/2 + u(j)'^/A, and K, r, u are parameters which can be determined in MD simulations by measuring the 
interface width ^ = \jK/r, the interfacial tension 7 = 2A/2r^^/3u, and the two homogeneous equilibrium phases 
= rkx/rju (= ±1 in our case). The two coupled equations of motion are the CH convection-diffusion equation for 



and the Navier-Stokes equation (with the addition of the capillary force density): 



'di 



(1) 



= -Vp V • (t" pV<t) + mpsext, (2) 

together with the incompressibility condition V • v = 0. Here M is the phcnomcnological mobility coefficient, pm is 
the fiuid mass density, p is the pressure, cr^ is the Newtonian viscous stress tensor, /iV0 is the capillary force density 
with n = 5F/S(p being the chemical potential, and mpgext is the external body force density (for Poiseuille flows). The 
boundary conditions at the solid surface are Vn = 0, = (n denotes the outward surface normal), a relax;ational 
equation for (j) at the solid surface: 

^+v^cP = -Tm, (3) 



- + (v.V)v 
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and the continuum form of the GNBC: 

pvfP = -r,dnV, + L{cj>)d,cj>, (4) 

Here F is a (positive) phenomenological parameter, I/(<^) = Kdncj) + d')wf{4>)/d(f) with 7u,/(<^) being the fluid-solid 
intcrfacial free energy density, /? is the slip coefficient, v^'^'p is the (tangential) slip velocity relative to the wall, r} is 
the viscosity, and L{(j))dx4> is the uncompensated Young stress. 

Compared to the Navier boundary condition, the additional component captured by the GNBC in Eq. (4) is the 
uncompensated Young stress aj^. Its differential expression is given by 

aj, = L{cf>)dxcj> = [KdncP + dj^f{cj>)/d4>]dxcP (5) 

at ^; = 0. Here the z coordinate is for the lower fluid-solid interface where dn = —dz (same below), with the under- 
standing that the same physics holds at the upper interface. It can be shown that the integral of this uncompensated 
Young stress along x across the fluid-fluid interface yields 

/ dx[L{cj>)dxcj>]{x,0)=jcose'/''f + Aj^f, (6) 

J int 

where J^^^ dx denotes the integration along x across the fluid-fluid interface, 7 is the fluid-fluid interfacial tension, 
^sur/ ^j^g dynamic contact angle at the solid surface, and A^^f is the change of ^wf{4>) across the fluid- fluid 
interface, i.e., A'jyjf = J ^^^dxdx'^w The Young's equation relates A^y^f to the static contact angle Of^'^^ at the 
solid surface: 

7coser'^ + A7^j = 0, (7) 

which is obtained as a phenomenological expression for the tangential force balance at the contact line in equilibrium. 
Substituting Eq. (7) into Eq. (6) yields 

/ dx[L{(l))dx(l)] {x, 0) = 7(cos el^""^ - cos (8) 

J int 

implying that the uncompensated Young stress arises from the deviation of the fluid-fluid interface from its static 
configuration. Equations (5), (6), (7), and (8) will be derived in Sec. V. In essence, our results show that in the 
vicinity of the MCL, the tangential viscous stress —rjdnVx as postulated by the usual Navier boundary condition can 
not account for the contact-line slip profile without taking into account the uncompensated Young stress. This is seen 
from both the MD data and the predictions of the continuum model. 

Besides the external conditions such as the shear speed V and the wall separation H, there are nine parameters in 
our continuum model, including p^, /?, 7: I0±|, M, F, and 9l"-''K The values of pm, 11, /?, ^, 7, |0±|, and 
were directly obtained from MD simulations. The values of the two phenomenological parameters M and F were fixc;(l 
by an optimized comparison with one MD flow fleld. The same set of parameters (corresponding to the same local 
properties in a series of MD simulations) has been used to produce velocity fields and fiuid-fluid interface shapes for 
comparison with the MD results obtained for different external conditions {V, H, and flow geometry). The overall 
agreement is excellent in all cases, demonstrating the validity of the GNBC and the hydrodynamic model. 

The CH hydrodynamic formulation of immiscible two-phase flow has been successfully used to construct a continuum 
model. We would like to emphasize that while the phase-field formulation does provide a convenient way of modeling 
that is familiar to us, it should not be conceived as the unique way. After all, what we need is to incorporate our key 
finding, the GNBC, into a continuum formulation of immiscible two-phase flow. The GNBC itself is simply a fact 
found in MD simulations, independent of any continuum formulation. 



IV. MOLECULAR DYNAMICS I 

A. Geometry and interactions 

MD simulations have been carried out for two-phase immiscible flows in Couette geometry (see Figs. 10 and 11) 
[23]. Two immiscible fluids were confined between two planar solid walls parallel to the xy plane, with the fluid-solid 
boundaries defined by ^ = 0, H. The Couette fiow was generated by moving the top and bottom walls at a constant 
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speed V in the ±.t directions, respectively. Periodic boundary conditions were imposed along the x and y directions. 
Interaction between fluid molecules separated by a distance r was modeled by a modified LJ potential 



Uff = 4e ((j/r)^^ - Sff {a/r 



where e is the energy scale, a is the range scale, with Sff = 1 for like molecules and Sff = —1 for molecules of different 
species. (The negative Sff was used to ensure immiscibility.) The average number density for the fluids was set at 
p = 0.81cr~^. The temperature was controlled at 2.8e/kB- (This high temperature was used to reduce the near-surface 
layering induced by the solid wall.) Each wall was constructed by two [001] planes of an fee lattice, with each wall 
molecule attached to a lattice site by a harmonic spring. The mass of the wall molecule was set equal to that of the 
fluid molecule m. The number density of the wall was set at pw = 1.860-"^. The wall-fluid interaction was modeled 
by another LJ potential 



C/w/ = 4e™/ {<^wf/r)^'^-Siufi(Ttnf/rf 



with the energy and range parameters given by e^/ = 1.16e and ayjf = 1.04(7, and Syjf for specifying the wetting 
property of the fluid. There is no locked layer of fluid molecules at the solid surface. We have considered two cases. 
In the symmetric case, the two fluids have the identical wall- fluid interactions with S^tf = 1. Consequently, the static 
contact angle is 90° and the fluid-fluid interface is flat, parallel to the yz plane. In the asymmetric case, the two fluids 
have different wall-fluid interactions, with S^f = 1 for one and S^f = 0.7 for the other. As a consequence, the static 
contact angle is 64° and the fluid-fluid interface is curved in the xz plane. In most of our simulations, the shearing 
speed V was on the order of O.lye/m, the sample dimension along y was 6.8a, the wall separation along z varied 
from H = 6.8a to 68a-, and the sample dimension along x was set to be long enough so that the uniform single-fluid 
shear flow was recovered far away from the MCL. Steady-state quantities were obtained from time average over 10'''t 
or longer where r is the atomic time scale •^/roa^/e. Throughout the remainder of this paper, all physical quantities 
are given in terms of the LJ reduced units (defined in terms of e, a, and m). 




FIG. 11. Schematic of simulation geometry, (a) Static configurations in the symmetric and asymmetric cases. Here fluid 2 
is sandwiched between fluid 1 due to the periodic boundary condition along the x direction. In the symmetric case, the static 
contact angle is 90° and the fluid-fluid interface is flat, parallel to the yz plane. In the asymmetric case, the static contact 
angle is not 90° and the fluid-fluid interface is curved in the xz plane, (b) Dynamic configuration in the symmetric case. 



B. Boundary- layer tangential wall force 

We denote the region within zq = 0.85a from the wall the boimdary layer (BL) (see Fig. 12). It must be thin 
enough to ensure sufficient precision for measuring the slip velocity at the solid surface, but also thick enough to 
fully account for the tangential wall-fiuid interaction force, which is of a finite range. The wall force can be singled 
out by separating the force on each fluid molecule into wall- fluid and fluid-fluid components. The fluid molecules in 
the BL, being close to the solid wall, can experience a strong periodic modulation in interaction with the wall. This 
lateral inhomogeneity is generally referred to as the "roughness" of the wall potential [3] . When coupled with kinetic 
collisions with the wall molecules, there arises a nonzero tangential wall force density that is sharply peaked at 
z ~ 2:0/2 and vanishes beyond z « (see Fig. 13). From the force density 5™, we define the tangential wall force per 
unit area as G^{x) = J^" dzg^{x, z), which is the total tangential wall force accumulated across the BL. 
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FIG. 12. Boundary layer at the lower fluid-solid interface. The empty and solid circles indicate the instantaneous molecular 
positions of the two fluids projected onto the xz plane. The solid squares denote the wall molecules. The dashed line indicates 
the level oi z = zq. 
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FIG. 13. Subdividing the BL into thin sections, we plot the reduced tangential wall force density as a function of distance z 
away from the boundary. The solid lines are averaged <?™(z) in thin sections at different x, normalized by the corresponding 
total wall force per unit area. The dashed line is a smooth Gaussian fit. It is seen that g^{z) is a function sharply peaked at 



zo/2. Note that dz[g^{x,z)/G^{x)] 



1. 



The short saturation range of the tangential wall force may be understood as follows. Out of the BL, each fluid 
molecule can interact with many wall molecules on a nearly equal basis. Thus the modulation amplitude (the 
roughness) of the wall potential would clearly decrease with increasing distance from the wall. That's why the 
tangential wall force tends to saturate &t z zq, which is still within the cutoff distance of wall- fluid interaction. 
On the contrary, the normal wall force arises from the direct wall-fluid interaction, independent of whether the wall 
potential is rough or not. Consequently, the normal wall force saturates much slower than the tangential component. 

We have measured the slip velocity and the tangential wall force in the BL. Spatial resolution along the x direction 
was achieved by evenly dividing the BL into bins, each Aa; = 0.85(7 or 0.425a". The slip velocity was obtained as 
the time average of fluid molecules' velocities inside the BL, measured relative to the moving wall; the tangential wall 
force was obtained from the time average of the total tangential wall force experienced by the fluid molecules in 
the BL, divided by the bin area in the xy plane. As reference quantities, we also measured G™*^ in the static {V ~ 0) 
conflguration. Figure 14 shows vf"^^ and measured in the dynamic configuration and G™° measured in the static 
configuration. It is seen that in the absence of hydrodynamic motion = 0), the static tangential wall force G^^ is 
not identically zero everywhere. Instead, it has some fine features in the contact-line region (a few cr's) (sec Fig. 14b). 
This nonzero static component in the tangential wall force arises from the microscopic organization of fiuid molecules 
in the contact-line region. 

The static component is also present in measured in the dynamic configuration, as shown by Fig. 14b. To see 
the effects arising purely from the hydrodynamic motion of the fliiids, we subtract G!^° from G™ through the relation 

' 

where G™ is the hydrodynamic part in G^. In the notations below, the over tilde will denote the difference between 
that quantity and its static part. We find the hydrodynamic tangential wall force per unit area, G^, is proportional 
to the local slip velocity w*'*^: 

G^{x) = -l3vfP{x), (9) 
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where the proportionahty constant is the shp coefficient. In Fig. 15, is plotted as a function of v^^'^p. The symbols 
represent the values of and t;^''^ measured in the BL. The lines represent the values of G^ calculated from — /Jv*'*^ 
using uj'*^' measured in the BL and (3 = {(3\p\ + fi^p^) / {pi +P2), with /3i^2 the slip coefRcients for the two fluid species 
and pi_2 the molecular densities of the two fluid species measured in the BL. Independent measurements determined 
Pi = P2 = l.2^/aji/a^ for the symmetric case, /3i = 1.2y^em/a^ and /?2 = for the asymmetric case. 

(The dependence of (3 on /3i,2 and pi^2 assumes the two fluids interact with the wall independently. The desired 
expression is obtained by expressing G^ as the weighted average of G^^ = — /3iw*''^^ and G!^^ = —(32Vx^^'^ noting 
^Mpi ^ ^sUpi within' 10%). 




FIG. 14. Slip velocity and tangential wall force (in reduced units) measured in the BL at the lower fluid-solid interface. 

(a) The slip velocity — + V \i, plotted as a function of x. The solid line denotes the dynamic symmetric case with 

V = 0.25^y e/m and H — 13.6a; the dashed line denotes the dynamic asymmetric case with V = 0.2-^ e/m and H = 13.6cr. 
The slip at the contact line {x ^ 0) is near-complete, i.e., v^''''' « V. (b) The tangential wall force is plotted as a function of x. 
The solid line denotes in the dynamic symmetric case; the dashed line denotes G™ in the dynamic asymmetric case. The 
dotted line denotes G^" in the static symmetric case; the dot-dashed line denotes G™° in the static asymmetric case. Note 
that G^° vanishes out of the contact-line region. 
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FIG. 15. 1/G™ plotted as a function of 1/wj''^ (in reduced units). Symbols are MD data measured in the BL at different x 
locations. The circles denote the symmetric case with V = 0.25-^ e/m and H = 13.6ct; the squares denote the asymmetric case 

with V = 0.2 s/ c/rn and H = 13.6ct. The solid and dashed lines are calculated from Eq. (9) for the symmetric and asymmetric 
cases, respectively, as described in the text. The lower-right data segment corresponds to the lower BL, whereas the upper-left 
segment corresponds to the upper BL. The slopes of the two dotted lines are given by /3i~2, which are proportional to the slip 
length. 



C. Tangential fluid force 

In either the static equilibrium state (where G™ = 0) or the dynamic steady state (where G^ / 0), local force 
balance necessarily requires the stress tangential to the fluid-solid interface to be the same on the two sides. Therefore, 
the hydrodynamic tangential fluid force per unit area, G^, must be proportional to the slip velocity wj''^*: 
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Giix)=Pvf^{x), (10) 

such that GI{x) + G^^{x) = according to Eq. (9). (The MD evidence for this force balance will be presented in Sec. 
VII.) Physically, is the hydrodynamic force along x exerted on a BL fluid element by the surrounding fluids, and 
may be expressed as 

Gl{x) = Jq° dz[dxaxx{x,z) + dzazx{x,z)] 
= ^zx{x, zo) + dx Jq" dzdxxix, z), 

using the fact that azx{x,0) = 0. (More strictly, azx{x,0~) = because there is no fluid below z = 0, hence no 
momentum transport across z = 0.) Here axx{zx) = (Txx(zx) — '^xx{zx)' ^^^^ being the normal component and ai^J 
the tangential component of the fluid stress tensor in the dynamic (static) configuration. 



D. Sharp boundary limit 

The form of in Eq. (11) is derived from the fact that the tangential wall force is distributed in a BL of finite 
thickness (Figs. 12 and 13). Now we take the sharp boundary limit by letting the tangential wall force strictly 
concentrate at z = 0: g^ix, z) = G'^{x)d{z) with the same G'^{x) per unit area. As shown in Fig. 13, the tangential 
wall force density is a sharply peaked function. By taking the sharp boundary limit the normalized peaked function 
is replaced by 6{z). Rewriting G^ in Eq. (11) as 

Gl{x) = /q- dz[dxd-xxix, z) + dzo-zxix, z)] 

= CFzxix, 0+) + /q+ dz[dxdxx{x, Z) + dzCFzxix, z)\, 

we obtain 

Gi{x)=azx{x,Q+)=Pvfnx), (12) 

because local force balance requires dx^xx + dzO'zx = above 2 = 0"*". Therefore, in the sharp boundary limit dzx 
varies from azx{x, 0~) = to ^zx{x, 0+) = Gl{x) at ^ = such that 

{V-a)-ii = Gl{x)6{z), 

in balance with the tangential wall force density g^{x,z) = G^{x)6{z). Equation (12) may serve as a boundary con- 
dition in hydrodynamic calculation if a continuum (differential) form of azx{x, 0+) is given. This will be accomplished 
in Sec. VB. 



V. CONTINUUM HYDRODYNAMIC MODEL 



For decades numerous models have been proposed to resolve the boundary condition problem for the contact- 
line motion [10-16], but so far none has proved successful by reproducing the slip velocity profiles observed in MD 
simulations [18,19]. In particular, based on the extreme (velocity) variations in the slip region, a breakdown of local 
hydrodynamic description in the immediate vicinity of the MCL has been suggested [19]. 

The main purpose of this paper is to present a continuum hydrodynamic model that is capable of reproducing MD 
results in the molccular-scalc vicinity of the MCL [23]. For this purpose, wc have derived a differential form for Eq. 
(12) (the continuum GNBC, Eq. (4)) using the Cahn-Hilliard hydrodynamic formulation of two-phase flow [15,16]. 
Our model consists of the convection-diffusion equation in the fluid-fluid interfacial region (Eq. (1)), the Navier-Stokes 
equation for momentum transport (Eq. (2)), the relaxational equation for the composition at the solid surface (Eq. 
(3)), and the GNBC (Eq. (4)). 



A. Cahn-Hillicird free energy functional 



The CH free energy was proposed to phenomenologically describe an interface between two coexisting phases [27]. 
In terms of the composition order parameter (p — {p2 — Pi)/{p2 +^1)7 the CH free energy functional reads 



-I 



dr 
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with /((/)) 



. Two thermally stable phases are given by (j)± = ±^/r/u at which df/d(j) = 0. An 



interface can be formed between the phases of and </)_ in coexistence. 



The chemical potential /x is defined by 



1. Chemical potential 



SF , 
/i = — = -KV (j) -r4> + u4 



from which the diffusion current J = — MV/x is obtained with M being the mobility coefficient. The convection- 
diffusion equation (Eq. (1)) comes from the continuity equation 

D(l) d<J) „ , 

2. Interfacial tension 

A few important physical quantities can be derived from the CH free energy. We first derive the interfacial tension 7 
for the interface formed between ^_|_ and (p- . In equilibrium, the spatial variation of (j) is determined by the condition 
that /Li(r) is constant, i.e., 

—Kd^(f) — r(f> + u(f>^ = constant. 

Here the interface is assumed to be in the xy plane with the interface normal along the z direction and the constant 
equals to zero because lim2_>±oo <t> = <t>± and lim^^±oo A* = 0- The interfacial profile is solved to be 



4>o{z) = <!>+ tanh 



with ^ = \fKjr being a characteristic length along the interface normal. The first integral is 

where the integral constant C equals f{4>±). It follows the interfacial free energy per unit area, i.e., the interfacial 
tension, is given by 



'y = j dz 

Using the interfacial profile (^0(2), we obtain 



-K{dAr+f{<t>)-f{cj^±) 



j dzK{d,4>f. 



1 



V2^ 



— J dz cosh ^ z 



3C 



3u 



3. Capillary force and Young stress 

We now turn to the forces arising from the interface. Consider a virtual displacement u(r) and the corresponding 
variation in (p, 5^{r) = — u(r) • \/(f). The change of the free energy due to this Sep is 



6F =Jdr 



df{ct>) 



+ Jdr 



d 



d{d, 



= Jdr [iiS<j)] +Jds [Kdr, 
= - Jdr[g-u] + Jds [crliUi] , 
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where g = /iV0 is the capillary force density in the Navier-Stokes equation (Eq. (2)), and a^^ = —Kdn4>di(j) is the 
tangential Young stress (the i direction is along the fluid-solid interface, i _L n). 

The body force g(r) = iiV4> can be reduced to the familiar curvature force in the sharp interface limit [28]. The 
unit vector normal to the level sets of constant (j) is given by m = V^/|V^| and 

/iV(^ = [-KW'^(j) -r(t> + u(t)^]\V(j)\m 

= -i4:V2^|V</)|m + i-Kdl^cj) -r<j) + u03]|v0|m 

where V* and dm denote the differentiations tangential and normal to the interface respectively. For gently curved 
interfaces, the order parameter (j) along the interface normal can be approximated by the one-dimensional stationary 
solution (^0, i-e., —Kd^<p — r(l) + ucj)^ w 0. Hence, /uV0 « —KWl<f)\W(l)\Toa., from which we obtain the desired relation 

where k = — (/)/|V(/'| is the curvature and 7 ~ / dlmK {V^))^ « / dlmK {dm'P)^ is the interfacial tension, with 1^ 
being the coordinate along the interface normal and the interface located at = 0. 

For gently curved interfaces, Kdn(j) ~ K dm(f> cos 9^"^^ ^ , where n is the outward (solid) surface normal, m the (fluid- 
fluid) interface normal, and O'^'^'^i the angle at which the interface intersects the solid surface (n • m = cos6'*"'"'^). For 
the tangential Young stress cr^ = Kdn4'9x(t> at 2; = where n = — z and i = x, the integral J.^^ dxaj^ along x across 

the interface equals to j^^^dxK^n4>^a:4> = {f^^^d(pKdm(t>) cos 9^'"''^ , where J.^^d^Kdm(f> = j^^^dlmK {dm4>f = 7- 
Hence, 

/ da;c7^^ =7cos6i™''-^, (13) 

J int 

where 9'^^^^ may be the dynamic contact angle at the solid surface 9^^^^ or the static contact angle 9^"^^ . This 
J.^^dxaJ^ is the tangential force per unit length at the contact line (aligned along y), exerted by the fluid- fluid 
interface of tension 7, which intersects the solid wall at the contact angle . So it equals to 7 cos . 



4- Young's equation 



The Young's equation for the static contact angle 0™'"/ can be derived as well. Consider the interfacial free energy 
at the fluid-solid interface, F^f = J ds'ywf{(j)). Minimizing the total free energy F + F^f with respect to (p at the 
solid surface yields 



from which an equation of local tangential force balance 



= 0, 



(^°x + ^xlwfifl^eq) =0, 



(14) 



(15) 



is obtained at 2; = 0. Here aj^ = aj,^ + dxjwf is the uncompensated Young stress (flrst introduced in Eq. (5)), (j)eq 
is the equilibrium composition field, and cr^^, denotes the static Young stress a^^{(l)eq). Integrating Eq. (15) along x 
across the interface leads to the Young's equation j cos9l'^'^-f + A7^/ = (Eq. (7)), where jcos9g^^^ = J^^^dxa°^ 
and A'y^jf = J-^^dxdx'ywf{<t>) is the change of fluid-solid interfacial free energy per unit area across the fluid- fluid 
interface. A microscopic picture for the Young's equation as an (integrated) equation of tangential force balance will 
be elaborated in Sec. VHB 1. 



B. Two boundary conditions 

Equations (14) and (15) are boundary conditions for the equilibrium state. In the dynamic steady state, however, 
neither Kdn(p + d-fwf{<j))/d(j) = L{(j)) nor cr^ + dx'jwf{(j)) = L{(j))dx<l> vanishes. In fact, the nonzero L{(f)) is responsible 
for the relaxation of <j) at the solid surface while the nonzero L{(l))dx4> is necessary to a slip boundary condition that 
is able to account for the near-complete slip at the MCL. 
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The convection-diffusion equation (Eq. (1)) is fourth-order in space. Consequently, besides the usual impermeability 
condition = 0, one more boundary condition is needed. The dynamics of (p at the solid surface is plausibly assumed 
to be relaxational, governed by the first-order extension of Eq. (14). More explicitly, when the system is driven away 
from the equilibrium, both d(j)/dt + v • and L{(j)) become nonzero, and they are related to each other by a linear 
relation 

|^+vV0cxL(0). 

This leads to Eq. (3) with T introduced as a phenomenological parameter. 
The GNBC (Eq. (4)) is obtained by substituting 

^.Ax, 0+) = a, Ax, 0) - (7°, (a;, 0) = a^Jx, 0) + aj^x, 0) - C7%{x, 0) 

= aUx,0) + <7l,{x,0) + d,j^f (16) 

into Eq. (12). Here the hydrodynamic tangential stress a^x is decomposed into a viscous component o"|,j, and a 
non- viscous component (t^. The viscous component is simply given by c!;'^ = r]dzVx] the non- viscous component 
is the uncompensated Young stress a^, given by = ct^ + dx'ywf{4') (Eq. (5)). According to Eq. (15), this 
uncompensated Young stress vanishes in the equilibrium state. But in a dynamic configuration, from the integral of 
along X across the fluid-fluid interface (Eqs. (6), (7), and (8)) 

/ dxu^x = 7 cos e^^^f + Aj^f = 7(cos ^^""^ - cos ef'^f), 

J int 

there is always a non- viscous contribution to the total tangential stress a^x as long as the fluid-fluid interface deviates 
from its static configuration. 

In Sec. VI we will show that the GNBC, with the uncompensated Young stress included, can account for the slip 
velocity profiles in the vicinity of the MCL, especially the near-complete slip at the contact line. In Sees. VII B and 
VII C we will present more MD evidence supporting the GNBC. A "derivation" of the GNBC, based on the tangential 
force balance (Sec. VII B) and the tangential stress decomposition (Sec. VII C), will be given in Sec. VIII. 



C. Dimensionless equations 



Dimensionless equations suitable for numerical computation are obtained as follows. We scale </> by \4>± 



length by ^ = \jK/r, velocity by the wall speed V, time by $,/V, and pressure/stress by "qV/S^. In dimensionless 
forms, the convection-diffusion equation is 

d(j) 



dt 



+ V • V(/) = £dV'(-V> + r ), 



the Navier-Stokes equation is 

n 

the relaxational equation for <^ at the solid surface is 

d(j) 

-g^ + Vxdx(p = -V. 

and the GNBC is 



dv 

^ + (v.V)v 



■ + <^')V.^, 



V2 
3 



Sn'/'- ^COs6'r''-^ S^(0) 



B 



%/2 



COS eT^ f sA<l>) 



d-r. 



dnVx- 



(17) 



(18) 



(19) 



(20) 



Here s^{(j)) = (7r/2) cos(7r^/2) is from the fiuid-soHd interfacial free energy 

7»/((/') = (A7^//2)sin(^0/2), 

which denotes a smooth interpolation between ±A7^j/2. Five dimensionless parameters appear in the above equa- 
tions. They are (1) Cd = Mr/V^, which is the ratio of a diffusion length Mr/V to (2) Tl = pV^jr], (3) 
B = r'^^/urjV = 3j/2\/2rjV , which is inversely proportional to the capillary number Ca = rjV/'-f, (4) V, = KT/V, and 
(5) £«((/)) =Tj/p{(j))^, which is the ratio of the slip length /s(0) =r]/f3{(p) to where /?(</)) = (1 - (^)/3i/2 (1 -F0)/32/2. 
A numerical algorithm based on a fixed uniform mesh has been presented in Ref. [23] . 
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VI. COMPARISON OF MD AND CONTINUUM RESULTS 



To demonstrate the validity of our continuum model, we have obtained numerical solutions that can be directly 

compared to the MD results for flow fleld and fluid-fluid interface shape. We have carried out the MD-continuum 
comparison in such a way that virtually no adjustable parameter is involved in the continuum calculations. This is 
achieved as follows. 

There arc totally nine material parameters in our continuum model. They are p„i, r], /3, ^, 7, \(f>±\, M, F, and Of^"^^ . 
(Note (1) For the asymmetric case, two unequal slip coefficients /3i and (32 are involved in /3; (2) The three parameters 
^, 7, and are equivalent to the three parameters K, r, and u in the CH free energy density; (3) is for 

^Jwf = — 7COs6'f"'"-^.) Among the nine parameters, seven are directly obtainable (measurable) in MD simulations. 
They are pm, Vi /3i,2, Ci 7: I0±l: ^-^d Of^^^ . (The fluid mass density Pm is set in MD simulations, the viscosity 77 
and the slip cocfhcients /3i^2 can be measured in suitable single-fluid MD simulations, the interfacial thickness ^ can 
be obtained by measuring the interfacial profile (/> = {p2 — Pi)/{p2 + Pi) in MD simulations, the interfacial tension 
7 can be obtained by measuring an integral of the pressure/stress anisotropy in the interfacial region [29], \4)±\ = 1 
means the total immiscibility of the two fluids, and the static contact angle Of^"^^ is directly measurable.) The two 
phenomenological parameters M and F have been introduced to describe the composition dynamics in the interfacial 
region. Their values are fixed by an optimized MD-continuum comparison. That is, one MD flow field is best matched 
by varying the continuum flow field with respect to the values of M and F. Once all the parameter values arc obtained 
(7 measured in MD simulations and 2 fixed by one MD-continuum comparison), our continuum hydrodynamic model 
can yield predictions that can be readily compared to the results from a series of MD simulations with different external 
conditions {V, H, and ffow geometry). The overall agreement is excellent in all cases, thus demonstrating the validity 
of the GNBC and the hydrodynamic model. We emphasize that the MD-continuum agreement has been achieved 
both in the molecular-scale vicinity of the contact line and far way from the contact line. This opens up the possibility 
of not only continuum simulations of nano- and microfluidics involving immiscible components, but also macroscopic 
immiscible flow calculations that are physically meaningful at the molecular level. (Molecular-scale details may be 
resolved through the iterative grid redistribution method without significantly compromising computation efficiency, 
see [30,31]). 

A. Immiscible Couette flow 

1. Two symmetric cases 

In Figs. 16 and 17 we show the MD and continuum velocity fields for two symmetric cases of immiscible Couette 
flow. In MD simulations, these two cases have the same local properties (fluid density, temperature, fluid-fluid 
interaction, wall-fluid interaction, etc) but different external conditions {H and V). Correspondingly, the continuum 
results are obtained using the same set of nine material parameters p^, f], f3 {= pi = /J2), ^, 7, |<?^±|) T, and 6f^^f . 

2. Two asymmetric cases 

In Figs. 18 and 19 we show the MD and continuum velocity fields for two asymmetric cases of immiscible Couette 
flow. In MD simulations, these two cases have the same local properties (fluid density, temperature, fluid-fluid 
interaction, wall-fluid interaction, etc) but different external conditions {H and V). Correspondingly, the continuum 
results are obtained using the same set of ten material parameters pm, P2, C) 7> \4>±\j Tj and 9f^^^ . In 

particular, among these parameters, [32 and Of^"^^ are measured in MD simulations while all the others directly come 
from the symmetric cases. Therefore, the comparison here is without adjustable parameters. 

3. From near- complete slip to uniform shear flow 

From Figs. 16, 17, 18, and 19, wc sec that at the MCL, the slip is near-complete, i.e., ~ and ~ V, 

while far away from the contact line, the flow field is not perturbed by the fluid-fluid interface and the single-fluid 
uniform shear flow is recovered. The slip amount in the uniform shear flow is 2lsV/{H + 2ls), vanishing in the limit 
oi H Ig. Here we encounter an intriguing question: In a mcsoscopic or macroscopic system, what is the slip profile 
which consistently interpolates between the near-complete slip at the MCL and the no-slip boundary condition that 
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must hold at regions far away? Large-scale MD and continuum simulations have been carried out to answer this 
question [31]. 

4- Steady-state fluid- fluid interface 

In Fig. 20 wc show the MD and continuum fluid-fluid interface profiles for one symmetric and one asymmetric cases 
whose velocity fields are shown in Figs. 16 and 18. 

0.05 I < < ■ • • < 1 




FIG. 16. Comparison of the MD (symbols) and continuum (lines) velocity profiles {v^ (x) and Vz (x) at different z levels) for 
a symmetric case of immiscible Couette flow {V = 0.25{e/m)^^^ and H = 13.6a). The profiles arc symmetric about the center 
plane z = H/2, hence only the lower half is shown aX z = 0.425cr (circles and solid lines), 2.125cr (squares and dashed lines), 
3.825(7 (diamonds and dotted line), and 5.525cr (triangles and dot-dashed lines). The MD velocity profiles were measured by 
dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85ct. 
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FIG. 17. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx{x) and Vz{x) at different z levels) for 
a symmetric case of immiscible Couette fiow {V = 0.05{e/m)''^'^ and H = 27.2a). The profiles are symmetric about the center 
plane z = H/2, hence only the lower half is shown at « = 0.85cr (circles and solid lines), 4.25cr (squares and dashed lines), 7.65cr 
(diamonds and dotted line), and 11.05cr (triangles and dot-dashed lines). The MD velocity profiles were measured by dividing 
the fluid space into 16 layers along z, each of thickness H/16 = 1.7a. 
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FIG. 18. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx[x) and ^^(a;) at different z levels) 
for an asymmetric case of immiscible Couette flow {V = 0.2(e/m)'^^^ and H = 13.6a), shown at z = 0.425cr (circles and solid 
lines), 2.975a (squares and loug-dashed lines), 5.525cr (diamonds and dotted line), 8.075(7 (up-triangles and dot-dashed lines), 
10.625cr (down-triangles and dashed lines), 13.175cr (left-triangles and solid lines). Although the solid lines are used to denote 
two different z levels, for each solid line, whether it should be compared to circles or left-triangles is self-evident (same for 
the next figure). The MD velocity profiles were measured by dividing the fluid space into 16 layers along z, each of thickness 
H/16 = 0.85a. 
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FIG. 19. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx (x) and Vz (x) at different z levels) for 
an asymmetric case of immiscible Couette fiow (V = 0.1{e/m)^^^ and H = 27.2a), shown at z = 0.85a (circles and solid lines), 
5.95a (squares and long-dashed lines), 11.05a (diamonds and dotted line), 16.15a (up-triangles and dot-dashed lines), 21.25a 
(down-triangles and dashed lines), 26.35a (left-triangles and solid lines). The MD velocity profiles were measured by dividing 
the fluid space into 16 layers along z, each of thickness H/16 = 1.7a. 
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FIG. 20. Comparison of the MD (symbols) and continuum (lines) fluid-fluid interface profiles, defined by pi = P2 {(p = 0). 
The circles and dotted line denote the symmetric immiscible Couette flow with V = 0.25(e/m)^''^ and H = 13.6a; the squares 
and dashed line denote the asymmetric immiscible Couette flow with V = 0.2(e/m)^''^ and H = 13.6f7. The MD proflles were 
measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85cr. 



B. Immiscible Poiseuille flow 



In order to further verify that the continuum model is local and the parameter values are local properties, hence 
applicable to different flow geometries, wc have carried out MD simulations and continuum calculations for immiscible 
Poiseuille flows. We find that the continuum model with the same set of parameters is capable of reproducing the 
MD results for velocity field and fluid-fluid interface profile, shown in Fig. 21. Similar to what we have observed in 
Couette flows, here at the MCL the slip is near-complete, i.e., Vx ^ and ~ V, while far away from the contact 

line, the flow fleld is not perturbed by the fluid-fluid interface and the single-fluid unidirectional Poiseuille flow is 
recovered. In particular, the slip amount in the unidirectional Poiseuille flow vanishes in the limit oi H ^ Ig. 

We emphasize that the overall agreement is excellent in all cases (from Fig. 16 to 21), therefore the validity of the 
GNBC and the hydrodynamic model is well affirmed. 




FIG. 21. Comparison of the MD (symbols) and continuum (lines) results for an asymmetric case of immiscible Poiseuille 
flow. An external force rngext ~ 0.05e/a is applied on each fluid molecule in the x direction, and the two walls, separated 
hy H = 13.6a, move at a constant speed V = 0.51(e/m)^^^ in the —x direction in order to maintain a time-independent 
steady-state interface, (a) Fluid-fluid interface proflles, defined by pi = p2 (0 = 0). (b) Vx{x) at different z levels. The profiles 
are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425cr (circles and solid line), 2.125a 
(squares and dashed line), 3.825a (diamonds and dotted line), and 5.525a (triangles and dot-dashed line). The MD profiles 
were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85a. 
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C. Flow in narrow channels 



It is generally believed that continuum hydro dynamic predictions tend to deviate more from the "true" MD results 
as the channel is further narrowed [32]. This tendency has indeed been observed but the deviation is not serious for 
H as small as 6.8a, as shown in Fig. 22. This deviation is presumably due to the short-range molecular layering 
induced by the rigid wall [33]. As the channel becomes narrower, the layered part of the fluids occupies a relatively 
larger space, thus making the MD-continuum comparison less satisfactory. 
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FIG. 22. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx{x) and Vz{x) at different z levels) for 
a symmetric case of immiscible Couette flow {V = 0.25(e/m)^/^ and H = 6.8(t). The profiles are symmetric about the center 
plane z = H/2, hence only the lower half is shown at z = 0.425(7 (circles and solid lines), f.275CT (squares and dashed lines), 
2.125cr (diamonds and dotted line), and 2.975(7 (triangles and dot-dashed lines). The MD velocity profiles were measured by 
dividing the fluid space into 8 layers along z, each of thickness H/8 = 0.85a. 



D. Temperature effects 



Most of our MD results have been obtained by setting the temperature at 2.8e/kB, above the liquid-gas coexistence 
region. Such a high temperature was used to reduce the fluid layering at the solid wall [33]. Similar two-fluid 
simulations have also been performed for temperatures ranging from 1.2e/fcB to 3.0e/fcB. We find that the MD results 
can always be reproduced by our continuum model, with material parameters (e.g. viscosity, interfacial tension, and 
slip length) varying with the temperature. In Fig. 23 we show the MD velocity profiles obtained at the temperatures 
1.4e/fcB and 2.8e/kB- It can be seen that they are qualitatively very close to each other. The quantitative difference 
is due to the different material parameters at different temperatures. 

Finally we list in Table I the parameter values in the continuum hydrodynamic model, used for the MD-continuum 
comparison at T = 2.8e/kB- 

TABLE I. Parameter values used in the continuum hydrodynamic calculations for the MD-continuum comparison at 
r = 2.8e/fes. 





Pm 


^ 0.8lm/a^ 


ri R 


i 1.95^im/(7^ 






Isl 


= rj/l3i « 1.3(7 


ls2 


= v/Pi ~ 1-3(7 or 3.3(7 








0.33(7 


7 R 


i 5.5e/a^ 


I0±I = 1 




M 


w 0.023(7*/Vrne 


F K 


i Q.QQa / ^Jmt 


cos6ir''^ =0 or 


« 0.38 
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FIG. 23. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx{x) at difTcrent z levels) for two 
symmetric cases of immiscible Couette flow at different temperatures, (a) The case of T = lAe/kB, V = 0.25{e/rny^'\ and 
H = 13.6(7, with a weak wall-fluid interaction and a high ratio of pn> to p (compared to case b here), (b) The case of T = 2.8e/fcB , 
V = 0.25{e/m)^^^ , and H = 13.6cr, with e-wf = 1.16e, a-wf = 1.04cr, 5wf = 1, and Pw/p = 2.3. The profiles are symmetric about 
the center plane z = H/2, hence only the lower half is shown aX z = 0.425cr (circles), 2.125<t (squares), 3.825cr (diamonds), and 
5.525cr (triangles). The MD velocity profiles were measured by dividing the fiuid space into 16 layers along z, each of thickness 
H/16 = 0.85a. The slip length for (a) is larger than that for (b). Therefore, far away from the contact line, the slip amount in 
(a) (w 0.1(e/m)^/^) is larger than that in (b) (w 0.05{e/my^^). 



VII. MOLECULAR DYNAMICS II 



We have formulated a continuum hydrodynamic model based on the CH free energy and the GNBC. The solutions 
of the model equations agree with the MD results remarkably well. This indicates that our model captures the right 
physics, and hence more MD evidences can be obtained to support the continuum GNBC (Eq. (4)) which takes into 
account the uncompensated Young stress. This necessarily requires a reliable measurement of fluid stress near the 
solid surface, plus a decomposition of the tangential stress into two components, one being viscous and the other 
interfacial, as expressed by Eq. (16). 



A. Measurement of fluid stress 



Irving and Kirkwood [34] have shown that in the hydrodynamic equation of momentum transport, the stress tensor 
(flux of momentum) may be expressed in terms of the molecular variables as 

tT(r, t) = (Tif (r, t) + a-u{r, t), 

where aK is the kinetic contribution to the stress tensor, given by 



mi 



rrii 



Y{v,t) 



^-V(r,f) 



TO,; 



<5(x,-r) , 



and <Tu is the contribution of intermolecular forces to the stress tensor, given by 



\ i j^i I 



Here m^, p^, and Xj are respectively the mass, momentum, and position of molecule i, V(r, i) is the local average 
velocity, Fj^ is the force on molecule i due to molecule j, and (• • •) means taking the average according to a normalized 
phase-space probability distribution function. 
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The Irving-Kirkwood expression has been widely used for stress measurement in MD simulations. However, as 
pointed out by the authors themselves [34], the above expression for au represents only the leading term in an 
asymptotic expansion, accurate when the interaction range is small compared to the range of hydrodynamic variation. 
As a consequence, this leading-order expression for cru is not accurate enough near a fluid-fluid or a fluid-solid interface. 
Unfortunately, this point has not been taken seriously. For the MCL problem, a knowledge of the stress distributions 
at both the fluid-fluid and the fluid-solid interfaces is of fundamental importance to a correct understanding of the 
underlying physical mechanism. Therefore, a reliable stress measurement method is imperative. 

To have spatial resolution along the x and z directions, the sampling region was evenly divided into bins, each 
A.T = 0.425(7 by Az = 0.85f7 in size. The stress components o", ,. and a^x were obtained from the time averages of the 
kinetic momentum transfer plus the fluid-fluid interaction forces across the fixed- a; and z bin surfaces. More precisely, 
we have directly measured the x component of the fluid- fluid interaction forces acting across the x{z) bin surfaces, in 
order to obtain the xx{zx) component of ctjj. For example, in measuring (tjjzx at a designated z-oricnted bin surface, 
we recorded all the pairs of fluid molecules interacting across that surface. Here "acting/interacting across" means 
that the line connecting a pair of molecules intersects the bin surface (the so-called Irving-Kirkwood convention [34]). 
For those pairs, we then computed auzx at the given bin surface from 

_ 1 

where Ssz is the area of z-oriented bin surface, (i, j) indicate all available pairs of fluid molecules interacting across 
the bin surface, with molecule i being "inside of zSsz" and molecule j being "outside of zSsz" (molecule i is below 
molecule j), and Fijx is the x component of the force on molecule i due to molecule j. For a schematic illustration 
see Fig. 24. 
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FIG. 24. Schematic illustration for the measurement of the zx component of cru- The horizontal solid lines (separated by 
short vertical lines) represent bin surfaces with surface normal along the z direction. Circles denote fluid molecules. The 
dashed lines connect pairs of interacting molecules. Here the bin surfaces and the molecules are projected onto the xz plane. 
Molecules that appear to be close to each other may not be in the interaction range if their distance along y is too large. A 
pair of interacting molecules may act across more than one bin surface. Here the (1,3) pair acts across the surfaces A and C 
while the (1,5) pair acts across the surfaces B and D. At each bin surface the stress measurement must run over all the pairs 
acting across that surface. For surface D, there are three pairs of interacting molecules (1,5), (2,4), and (2,5) that contribute 
to the zx component of au- 

B. Boundary-layer tangential force balance 

From the data of stress measurement, we now present the MD evidence for the BL tangential force balance, flrst 
introduced in Sec. IV C for obtaining Eq. (10) from Eq. (9). 

1. Static tangential force balance 

We start from the tangential force balance in the static configuration {V = 0). As first pointed out in Sec. IV B, 
the static tangential wall force G™" shows molecular-scale features in the contact-line region, due to the microscopic 
organization of fluid molecules there. Then according to the local force balance, static fluid stress must vary in such 
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a way that the total force density vanishes. An integrated form of the static tangential force balance is given by 
G^°ix) + Gl°{x) = 0, where G^"(a;) = J^" dzg}^°{x,z) and the static tangential fluid force Gl°{x) is of the form 

Gl'^ix) = / dz[d,a°,{x, z) + d^aUx, z)] = a^^, zo) + / dzal,{x, z). (21) 
JO Jo 

Here and are the xx and zx components of fluid stress in the static configuration, both measured as reference 
quantities. In Fig. 25 we show dza'^j.^, a^^{zo), Gl^, and G™" (which is the same as in Fig. 14b). In the symmetric 
case, fi.T(j"j.(x, Zo) J.^^dxdx Jq° dza^^{x,z), and J^^fdx dzgJl^^{x, z) all vanish because = 90°. For the 

asymmetric case, G™°(a;) + G|°(x) = means 

a%{x, zo) + / dza°,{x, z) + / dzg^\x, z) = 0, (22) 
Jo Jo 

which corresponds to the continuum equation (Eq. (15)) \aj^ + dxjwA = at the solid surface. Here a'^JxjZo) 
corresponds to the continuum <jY-,.{'i'eq) at the solid surface while dx /q " dza^^+ " dzg'^^ corresponds to the continuum 
dxlmfifpeq) at the solid surface. The Young's equation (Eq. (7)) 7 cos 6'^"''-'' + = is then obtained through 

integration, using 

/ dxal, {x,zo)= [ dxal, = 7 cos 9T'f (23) 

J int J int 



and 

r-Zo 



/ dxdx / dzal^ + dx dzg^° = / dxdx^v^f = Aj^f. (24) 

Jint Jo J int Jo J int 

Here 7 cos Of^''^ and A'y^f arc the two tangential forces per unit length along the contact line (along y) , the former 
due to the tilt of the fluid- fluid interface (6'|"'"'^ ^ 90°) while the latter due to the different wall- fluid interactions for 
the two fluid species. In fact, equations (23) and (24) are the microscopic definitions for the two continuum quantities 
"fcosOf^^^ and A^fwf in the Young's equation, whose validity is based on the microscopic tangential force balance 
expressed in Eq. (22). 




FIG. 25. Profiles of Jj'" dza'ix, <^'ix{zo), , and for the lower BL. The horizontal axes are x/a. We show dza^^ 
(e/cr^) in (a) for the symmetric case and in (c) for the asymmetric case. For clarity, a^x has been vertically displaced such that 
in the symmetric case, — far from the interface, and in the asymmetric case, a'^^ = at the center of the interface. The 
profiles of cr°j.(2;o), G^", and G™*^ (e/o"^) are plotted in (b) for the symmetric case and in (d) for the asymmetric case. The squares 
denote (t°x{zo), the diamonds denote Gl", and the solid triangles denote G^°. Both (b) and (d) show G^°(x) = —Gl°{x). 
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2. Dynamic tangential force balance 



An integrated form of the dynamic tangential force balance is given by G^'(a;) + Gl{x) = 0, where G^^{x) = 
Jq° dzg'^{x, z) and the dynamic tangential fluid force Gl{x) is of the form 

GI{x) = dz[da:(7xx{x,z)+dzaza:{x,z)]=a2a:{x,ZQ)+dx dzaxx{x,z). (25) 
Jo Jo 

Here and a^x axe the xx and zx components of fluid stress in the dynamic configuration. In Fig. 26 wc 
show Jq° dzaxx, o'zxizo), and G™. From a comparison between Figs. 25 and 26, we can see that the dynamic 
quantities dzcFxx, (^zx{zq), and indeed show features seen from the static quantities j^^" dza''^^, (7"^(zo), and G^°, 
respectively. This is particularly evident in the asymmetric case where the static quantities vary more appreciably. 
The reason to take the static quantities as reference quantities is now clear: a hydrodynamic quantity must be obtained 
from the corresponding dynamic quantity by subtracting its static part, formally expressed as 

Q dynamic l-^^ static ' 

where the over tilde denotes the hydrodynamic quantity (first introduced for G^ in Sec. IV B). In Fig. 27 we show 
the MD evidence for the BL hydrodynamic tangential force balance, which is expressed as G!^(a;) + Gl{x) = 0. This 
equation is necessary for Eqs. (9) and (10) to hold simultaneously. 

In summary, to verify the static/dynamic tangential force balance, we need to (1) identify the BL where the 

tangential wall force G™^"^ is distributed; (2) measure the normal and tangential components of stress cJxJ and CzJ 
according to the original definition of stress; (3) calculate the tangential fluid force gI^^^ according to Eq. (21)/(25). 




FIG. 26. Profiles of dza.. 
(a) for the symmetric case {V = 
The profiles of azx (20) and 



-20 -10 10 20 

•,x, cTzxizo), and for the lower BL. The horizontal axes are x/a. We show " dza^x {e/a^) in 

= 0.25-^ e/m and H = 13. 60-) and in (c) for the asymmetric case {V = 0.2y/ e/m and H = 13.6cr). 
{e/cr^) are plotted in (b) for the symmetric case and in (d) for the asymmetric case. The squares 



denote azx{zo) and the diamonds denote G™. Note that in the vicinity of the contact line, G^ + cJzx(zo) ^ 0. The importance 
of the a-gradient of the 2:-integrated normal stress dx Jq° dzaxx is therefore evident. 
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FIG. 27. Hydrodyriamic force balance for the lower BL. The circles represent the symmetric case {V — 0.25y^ e/m and 
H — 13.6a); the squares represent the asymmetric case (V = 0.2-^ e/m and H = 13.6a). The empty symbols denote G™; the 
solid symbols denote — G^. It is seen that G^{x) = —G{.{x). 

C. Tangential Young stress 



Here we present the MD evidence for the decomposition of the tangential stress. In a dynamic configuration, away 
from the interfacial region the tangential viscous stress cr"^ = r]{dzVx + dxVz) is the only component in the (single-fluid) 
tangential stress azx- But in the (two- fluid) interfacial region, the tangential stress azx can be decomposed into a 
viscous component a^^ and a non- viscous component cr^: 



V I 



Y 
zx ' 



where a^^ is still r]{dzVx + dxVz) and cr^ is the tangential Young stress, satisfying 

dx(7j,j.{x^ z) = ^ COS 9 d{z). 



(26) 



(27) 



Here 0d{z) is the dynamic interfacial angle at level z. Equation (26) is essential to obtaining Eq. (16) for dzx{x,0). 
In a static configuration, the viscous stress a^^ vanishes and a^x becomes cr°^, satisfying 



Eg = / dxa1^{x,z) = J cos ds{z), 

J int 



(28) 



where 6s{z) is the static interfacial angle at level z. Figure 28 shows that both (t^ and cr^^ are nonzero in the 
interfacial region only. The inset to Fig. 28 shows the evidence for Eqs. (27) and (28), which identify ct^ and a^x as 
the dynamic and static Young stresses. 




FIG. 28. Dynamic and static Young stresses at z = zo. The solid circles denote a'^^ in the static symmetric case; the empty 
circles denote ct^ in the dynamic symmetric case {V = 0.25-^ e/m and H = 13.6a); the solid squares denote a'^^ in the static 
asymmetric case; the empty squares denote cr^ in the dynamic asymmetric case {V = 0.2^ e/m and H = 13.6<t). Here aYx 
was obtained by subtracting the viscous component ri{dzVx + d-xV-S) from the total tangential stress a^^- Inset: Sd,s plotted as 
a function of ')cosdd,s at different z levels. Here 9d,s was measured from the time-averaged interfacial profiles. The symbols 
have the same correspondence as in the main figure. The data indicate T,d,s = 7cos^d,3. 
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Equations (27) and (28) can be derived from the mechanical definition for the interfacial tension 7 [29] : 



7 



= / dim [P±{lm) - P\\{1. 
J int 



i.e., 7 is the integral (along the interface normal across the interface) of the difference between the normal and 
parallel components of the pressure, where Im is the coordinate along the interface normal m, and P± and Py are the 

pressure-tensor components normal and parallel to the interface, respectively. (Note that far away from the interface 
the pressure is isotropic and P± = P\\)- When the interfacial angle 64 (or Og) is 90°, the interface normal m = x and 
the non- viscous stress tensor in the interfacial region is diagonal in the xyz coordinate system: 



' non— VISCOUS 



-P± 
-P|| 








-Pll 



-P±I+(P±-P||)(I-mm), 



where I is the identity matrix. According to this expression, when the interfacial angle 6d (or 6s) deviates from 
90° (see Fig. 29), the Young stress cr^ (or a^^) arises from the interfacial stress anisotropy as the off-diagonal zx 
component of the microscopic stress tensor: 



Z • (Tnon-viscous • X = (Pj_ - P|| ) [z • (I - mm) • x] 

-(P± - P\\)TnzTnx = (Pj_ - P||) cos^d(5) sin0d(5). 



where m. 



cos 6d(^s) and = sin O^^s) ■ It follows that 

.no), 



= lint ^^^^^ ^) = lint dx{P± - P|| ) COS 6»d(s) sin ^d(^) 

= fint '^lni{P± - P||) C0s6»d(^) 

= 7C0S 9d(s), 



where 0d{s) is treated as a constant along x and dx sin 9d(^s) = dim- 
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FIG. 29. Schematic illustration for the origin of the tangential Young stress cr^"^ as an off-diagonal component of the 
microscopic stress tensor. 



VIII. THE GENERALIZED NAVIER BOUNDARY CONDITION 



We want to "derive" the continuum GNBC using the MD results in Sees. IV and VII. For this purpose we first need 
to establish the correspondence between the stress components measured in MD and those defined in the continuum 
hydrodynamics. This correspondence is essential to obtaining the microscopic dynamic contact angle O^d^^ , which 
is defined in the continuum hydrodynamics (see Eq. (6) and (8)) but not directly measurable in MD simulations 
(because of the diffuse BL). 



A. MD-continuum correspondence 



It has been verified that for a BL of finite thickness, the GNBC is given by 

pvf^ix) =Gi(x)-Gi%x) 

^ dx Jo C'^^^^^' ^) ~ Z)] + [(^zx{x, Zo) - (T°a:{x, ^o)] , ^^^^ 
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in which only MD measurable quantities are involved. Now we interpret these MD-nicasured quantities in terms of 
the various continuum variables in the hydrodynamic model. In so doing it is essential to note the following. 

(1) a^x can be decomposed into a molecular component and a hydrodynamic component: axx = Txx + '^xx'- Mean- 
while, cr°^ can be composed into the same molecular component and a hydrostatic component: 

Physically, Txx is the normal stress (t°^ for the case of a flat, static fluid-fluid interface. Such an interface exists in 
the symmetric case {Of^^^ = 90°) for any value of H. It also exists in the asymmetric case [Of^"^^ ^ 90°) for i7 ^ oo 
(with vanishing curvature ~ 1/iJ). In either case, the interface has zero curvature and the hydrostatic stress a^x 
vanishes according to the Laplace's equation: cr°^ Txx as cTxx ~^ 0- (To be more precise, jcos9g"'^^ and Aj^f 
should be defined in Eqs. (23) and (24) when the fluid-fluid interface has zero curvature. However, in the asymmetric 
case where Ajyjf (= —jcosO^'^'^-f) is nonzero, there is a static interfacial curvature ~ 2cos9l^^^ /H. This results in 
a hydrostatic cr^x ^ which should be subtracted from a^x in the left-hand side of Eq. (24). That is, A-y^jf should 
be obtained with crj!^ replaced by Txx- Meanwhile, due to the static interfacial curvature, the static interfacial angle 
9a{zo) determined by cos6s{zo) = J^^^. dxa^^i^^ ^o) is a bit different from the true Of^^^ . In fact, cos9s{zo) deviates 
from cos^™''-^ by the BL- integrated curvature w 2 zq cos 91^'^ ^ /H.) The molecular component Txx exists even if there 
is no hydrodynamic fluid motion or fluid-fluid interfacial curvature. On the contrary, the hydrodynamic component 
arises from the hydrodynamic fluid motion and interfacial curvature. In the static configuration, cr^^ becomes 
cr^x ^ which comes from the interfacial curvature. 

(2) azx{x,ZQ) can be decomposed into a viscous component plus a Young component : (Jzx{x-iZq) = tT^j.(a;, 2:0) + 
a^xix, zq) with cr^^ = 7]{dzVx + dxV^) and /^^^ dxa^xix, Zq) = 7 cos 6*^(20) (Eqs. (26) and (27)). 

(3) cr"^(.T, Zq) is the static Young stress: i.e., J.^^^ dxa^xix, zo) = jcos9s{zo) (Eq. (28)). 

Using the above relations, we integrate Eq. (29) along x across the fluid-fluid interface and obtain 



Lt dxpvfP{x) = A [C dza^x"" ix, z)] + dxa^^x (x, zo) + 7 cos ^^(^o) 



where A 



-A [/;» dzag^{x,z)] -jcos9s{zo), 
Jq° dzaxJ^^^^^ is the change of the 2;-integrated Uxx"^^^^ across the interface: 



(30) 



A 



/ dza^^^''^^ = / dxdx \ dza'^^^''^\ 

Jo i Jint Jo 

According to the Laplace's equation, the hydrostatic stress is directly related to the static curvature Kg' 

-^^xx = 7K., 

and the ^-integrated curvature JJ^" dzns equals to cos ^5(20) — cos^™''-^. Hence, 

-A / dza^x^ = 7 / dzK, = 7 [cos 61^(20) - cos6»™''^] . (31) 
Jo Jo 

Substituting Eq. (31) into Eq. (30) yields 

/ dxpvfP {x) = A f " dza^x" {x, z) + f dxa^x {x, 20) + 7 cos 9d{zo) - 7 cos ^f''-^ . (32) 

J int Jo J int 

In order to interpret Eq. (32) in the continuum hydrodynamic formulation with a sharp BL, it is essential to note the 
following. 

(1) The sum of the first three terms on the right-hand side of Eq. (32) is the net fluid force along x exerted on the 
three fluid sides of a BL fluid element in the interfacial region. 

(2) The last term in the right-hand side of Eq. (32), —7 cos is the net wall force along x, A-f^^j = j^^^^ dxdx'^/wf, 
which arises from the wall-fluid interfacial free energy jump across the fluid-fluid interface, in accordance with the 
Young's equation A7i„/ -|- 7 cos ^J'"'-^ = 0. 

B. Extrapolated dynamic contact angle 

Now we take the sharp boundary limit to relate the net fluid force 



A f dzaxf {x, z)+ [ dxalx {x, 2^0) + 7 cos ^^(^o) 

Jo J int 
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in Eq. (32) to the tangential stresses (viscous and non-viscous) at the sohd surface. The purpose of doing so is 
to obtain the surface contact angle through extrapolation. Note that O^J^^^ is not directly measurable in MD 

simulations due to the diffuse BL. Only the extrapolated O^J^^^ can be compared to the contact angle in continuum 

calculations. 

In Sec. IV D, we take the sharp boundary limit by assuming a tangential wall force concentrated at z = 0: 
g'^{x,z) = G'^{x)S{z). While 5^ becomes a 5 function, G'^{x) per unit area remains the same. Using the equation 

of local force balance dxC^xx + dz<yzx — above z = 0^, we obtain (t^2;(.t,0^) = Gl{x) as the tangential stress at the 
solid surface (Eq. (12)). The extrapolation here follows this spirit. We turn to the Stokes equation in the BL: 

-dxP + dxcr^x + dzcr'zx + l^dx(t> = 0, (33) 

obtained from the a;-component of Eq. (2) by dropping the inertial term and the external force term. Integrating Eq. 
(33) along z across the BL and then along x across the fluid-fluid interface, we obtain 



^{/o"" dz[-p{x,z)+al^{x,z)]] + j^^^dxal^{x,ZQ) cose d{zo) 
= J^^^dxaUx,O)+^cos0l-^f. 

Two relations have been used in obtaining Eq. (34): 

(1) The capillary force density in the sharp interface limit [28] is given by 

^lda;(f) ~ Jk6{x - Xint), 

where k is the interfacial curvature and Xint the location of the interface along x (see Sec. V A3). 

(2) The ^;-integrated curvature gives 



(34) 



/" 

Jo 



dzK = cos9d{zo) — cos9^ 



surf 



The local force balance along x is expressed by Eq. (33). Accordingly, the tangential force balance for the BL fluids 
in the integration region (/.^^ dx dz) is expressed by Eq. (34), where A { j^" dz [—p{x, z) + cr^^i^, z)]} is the net 
fluid force on the left and right (^x-oriented) surfaces, /^^^ dxa^^i^' ^0) + 7 cos 9d{zo) is the fluid force on the z = zq 
surface, and J^^^ dxa^^{x,0) + jcose^^^ is the tangential fluid force at the solid surface. 
Substituting Eq. (34) into Eq. (32) and identifying the normal stress —p + a^^ with c^''-': 



we obtain 



/ dxpvfP (x) = / dxal^ (a;, 0) + 7 cos 9'J'''^ - 7 cos (9^'^^ , (35) 

J mt J int 

which is identical to the integral of the continuum GNBC along x across the fluid-fluid interface (Eqs. (4), (5), (6), 
and (8)). By doing so, we have assumed the tangential wall force is concentrated at z = 0. (This leads to Eqs. (33) 
and (34) between z = 0+ and z = zq.) In essence, the above extrapolation is to obtain the tangential stresses (viscous 
and non-viscous) at the solid surface when the limit of g^ixjz) — {x)d{z) is taken, as illustrated in Fig. 30. For 
G'^{x) distributed in the diffuse BL, there is actually no tangential stress at the solid surface. Only in the sharp 
boundary limit does a nonzero tangential stress appear at 2; = 0+, equal to the net fluid force accumulated from 
z = 0~ to z = zo the diffuse BL. 



(a) 



Z=Zo 



-z=0 



(b) 

nonzero 
tangential stress 



B 



zero ^ 

tangential stress 

wall force distributed inside 



D 



-z=0 



wall force concentrated at z=0 

FIG. 30. (a) For the real tangential wall force continuously distributed between z = and z = zq, the tangential stress 
in the fluid is continuous, vanishing at z = 0. (b) In the sharp boundary limit, the tangential wall force is considered to be 
concentrated at 2 = 0. Accordingly, the net fluid force on the fluids bound by A, B, C, and D (at « = 0"*") is considered to be 
zero. In other words, the net fluid force on the three surfaces A, B, and C is fully transmitted to the tangential stress at the 
surface D at 2; = 0"*" (Eq. (34)). So there arises an abrupt change of the tangential stress from at 2; = 0~ to at 2 = 0"*", 
by which the concentrated wall force is balanced. 
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It is worth emphasizing that the right-hand side of Eq. (35) is J^^^ dxa^xix, 0+), where a-^xix, O"*") is the extrapolated 
tangential stress in Eq. (12). Its continuum expression is given by Eq. (16). Equation (35) concludes our "derivation" 
of (an integrated form of) the continuum GNBC from the MD results presented in Sees. IV and VII. 

Equation (34) constitutes the basis for obtaining the dynamic contact angle O'^^^^^^ from MD data. The dominant 
behavior of ° dz {—p + a^^) = J^" dza^J^ is a sharp drop across the fluid-fluid interface. This stress drop implies a 
large curvature in the BL, which pulls the extrapolated 0^^^ closer to 0™''-^. We show the BL-integrated normal stress 
dzcfxx = Jq° dz [axx — o'xx] Fig. 31, where the large stress drop across the fluid- fluid interface is clearly seen. 
(In the asymmetric case, due to the small difference between and Txx, ^xx is not precisely the hydrodynamic stress 
a^J^. In fact, axx = (^xJ^ — (^xx ■ Nevertheless, A dza^£ is on the order of 2'yzocos6g'^'^^ /H, much smaller than 
the magnitude of the stress drop shown in Fig. 31.) In the partial-slip region at the vicinity of the MCL, J^" dzdxx 
shows a fast variation along x as well. This means that the BL tangential force balance cannot be established unless 
the gradient of the BL-integrated normal stress is taken into account (see Eqs. (11), (21), and (25)). It is worth 
pointing out that the stress variation depicted in Fig. 31 is also produced by the continuum hydrodynamic calculations, 
in semi-quantitative agreement with the MD data. In Fig. 32 we show the continuum profiles of layer-integrated 
pressure, obtained for a symmetric case of Couette flow. It is readily seen that the magnitude of the pressure change 
across the fluid-fluid interface decays away from the solid surface quickly. This conforms to the interface profile shown 
in Fig. 20: the interfacial curvature quickly decreases with the increasing distance from the solid surface. 
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FIG. 31. Jg^" dzaxx{x,z) = dz [axx{x,z) — a''^j.{x,z)j plotted as a function of x. The circles denote the symmetric case 

{V = 0.25-\/ e/m and H = 13.6cr) and the squares denote the asymmetric case {V = '^fl^fefm and H = 13.6cr). For clarity, 
was vertically displaced such that a^^ = far from the interface in the symmetric case, and for the asymmetric case, (t°^ = 
at the center of the interface (same as in Fig. 25). 
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FIG. 32. J^^^°l2 dzp{x, z) obtained for the symmetric case of V = 0.25^ 



-^o/2 



: 13.6cr, plotted as a function of x. 

The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425c7 (solid line), 
2.125(7 (dashed line), 3.825(7 (dotted line), and 5.525(7 (dot-dashed lines). The solid line denotes the BL-integrated pressure, 
in scmi-quautitativc agreement with the dzaxx{x, z) profile (circles) shown in Fig. 31. (Note that in Fig. 31 the MGL is 
placed at X « while here it is shifted to a: « —3.9(7, same as in Figs. 16, 20, and 23b.) 
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C. The importance of the uncompensated Young stress 

According to Eq. (26), the hydro dynamic tangential stress azx{x, zq) can be decomposed into the viscous component 
al^ and the non- viscous component aj^: 

^zx{x, zo) = a-l^{z) + arl^{x, Zq). 

In Fig. 33 we show that away from the interfacial region the tangential viscous stress zq) = 'q{dzVx-\-dxVz){x, zq) 

is the only nonzero component, but in the interfacial region ct^ = cjzx — <^zx ~ '^zx — '^Yx ~ '^zx is dominant, thereby 
accounting for the failure of the Navier boundary condition to describe the contact line motion. Therefore away from 
the MCL region the Navier boundary condition is valid, but in the interfacial region it clearly fails to describe the 
contact line motion. 

We have measured the z-integrated axx = <7xx — o'xx ™ the BL. The dominant behavior is a sharp drop across the 

interface, as shown in Fig. 31 for both the symmetric and asymmetric cases. As already pointed out in Sec. VIII B, 
this stress drop means a large curvature in the BL, which pulls the extrapolated 0^^' ^ closer to O^^^^ . The value of 
^T^^ obtained through extrapolation is 88 ± 0.5° for the symmetric case and 63 ± 0.5° for the asymmetric case at the 
lower boundary, and 64.5 ± 0.5° at the upper boundary. These values are noted to be very close to Of'^'f . Yet the 
small difference between the dynamic and static contact angles is essential in accounting for the near-complete slip at 
the MCL. 

In essence, om results show that in the vicinity of the MCL. the tangential viscous stress (j!^^. as postulated by 
the Navier boundary condition can not give rise to the near-complete MCL slip without taking into account the 
tangential Young stress in combination with the gradient of the (BL- integrated) normal stress <Jxx- For the static 
configuration, the Young stress is balanced by the normal stress gradient, leading to the Young's equation. It is only 
for a MCL that there is a component of the Young stress which is no longer balanced by the normal stress gradient, 
and this uncompensated Young stress is precisely the additional component captured by the GNBC but missed by 
the traditional Navier boundary condition. 
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FIG. 33. Two components of the hydrodynamic tangential stress aX z = zo, plotted as a function of x. The symbols connected 
by dashed linos denote aj^; the solid and dotted lines represent the viscous component. Here the symmetric case is represented 
by circles and solid line; the asymmetric case represented by squares and dotted line. In the contact line region the non- viscous 
component is almost one order of magnitude larger than the viscous component. The difference between the two components, 
however, diminishes towards the boundary, z = 0, due to the large interfacial pressure drop (implying a large curvature) in the 
BL, thereby pulling 9'/''^ closer to er""^' . 



IX. CONCLUDING REMARKS 



In summary, we have found the slip boundary condition, i.e., the GNBC, for the contact-line motion. Based on 
this finding, we have formulated a hydrodynamic model that is capable of reproducing the MD results of slip profile, 
including the near-complete slip at the MCL. It should be noted, however, that the present continuum formulation 
can not calculate fluctuation effects that are important in MD simulations. 

Based on the results and experiences obtained in the present study of the MCL, we have been working on the 
following. 

(1) Large scale MD simulations on two-phase immiscible flows show that associated with the MCL, there is a very 
large 1 /x partial-slip region, where x denotes the distance from the contact line. This power-law partial-slip region has 
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been reproduced in large-scale adaptive continuum calculations [30,31] based on the local, continuum hydrodynamic 
formulation presented in this paper. MD simulations and continuum solutions both indicate the existence of a universal 
slip profile in the Stokes- flow regime, well described by v'''''^^ (x) /V = 1/(1 + x/alg), where w^'^^ is the slip velocity, 
V the speed of moving wall, Ig the slip length, and a is a numerical constant ~ 1 [31]. A large 1/x partial slip 
region is significant, because the outer cutoff length scale directly determines the integrated effects, such as the total 
steady-state dissipation. While in the past the 1/x stress variation away from the MCL has been known [7], to our 
knowledge the fact that the partial slip also exhibits the same spatial dependence has not been previously seen, even 
though the validity of the Navier boundary condition at high shear stress has been verified [2-5] . 

(2) By explicitly taking into account the long-range wall-fluid interactions, our hydrodynamic model for two-phase 
immiscible flow at the solid surface can be used to investigate the dry spreading of a pure, nonvolatile liquid, attracted 
towards the solid by long-range van der Waals forces [9]. The precursor film, driven by the gradient of disjoining 
pressure due to the long-range force [35], was observed decades ago [36-38]. Nevertheless, the theoretical analysis 
remains to be difficult, because of a wide separation of different length scales involved [9,39-41]. Application of our 
model to this multi-scale problem is currently underway. 

(3) Decades ago it became well known that the driven cavity flow is incompatible with the no-slip boundary 
condition, since the latter would lead to stress singularity and infinite dissipation (known as corner-fiow singularity) 
[1,42]. While MD studies [43] have clearly demonstrated relative fluid- wall slipping near the corner intersection, the 
exact rule that governs this relative slip has been left unresolved. Based on the simulation technique developed in 
the present study, we have verified the validity of the Navier boundary condition in governing the fluid slipping in 
driven cavity flows. We have used this discovery to formulate a continuum hydrodynamics whose predictions are in 
remarkable quantitative agreement with the MD simulation results at the molecular level [26]. 
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